Three-dimensional relativistic simulations of rotating neutron-star collapse to a Kerr black hole 
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We present a new three-dimensional fully general-relativistic hydrodynamics code using high-resolution 
shock-capturing techniques and a conformal traceless formulation of the Einstein equations. Besides presenting 
a thorough set of tests which the code passes with very high accuracy, we discuss its application to the study 
of the gravitational collapse of uniformly rotating neutron stars to Kerr black holes. The initial stellar models 
are modelled as relativistic polytropes which are either secularly or dynamically unstable and with angular ve- 
locities which range from slow rotation to the mass-shedding limit. We investigate the gravitational collapse 
by carefully studying not only the dynamics of the matter, but also that of the trapped surfaces, i.e. of both 
the apparent and event horizons formed during the collapse. The use of these surfaces, together with the dy- 
namical horizon framework, allows for a precise measurement of the black-hole mass and spin. The ability to 
successfully perform these simulations for sufficiently long times relies on excising a region of the computa- 
tional domain which includes the singularity and is within the apparent horizon. The dynamics of the collapsing 
matter is strongly influenced by the initial amount of angular momentum in the progenitor star and, for initial 
models with sufficiently high angular velocities, the collapse can lead to the formation of an unstable disc in 
differential rotation. All of the simulations performed with uniformly rotating initial data and a polytropic or 
ideal-fluid equation of state show no evidence of shocks or of the presence of matter on stable orbits outside the 
black hole. 

PACS numbers: 04.25.Dm, 04.40.Dg, 04.70.Bw, 95.30.Lz, 97.60.Jd 



I. INTRODUCTION 

The numerical investigation of gravitational collapse of ro- 
tating stellar configurations leading to black-hole formation is 
a long standing problem in numerical relativity. However, it 
is through numerical simulations in general relativity that one 
can hope to improve our knowledge of fundamental aspects 
of Einstein's theory such as the cosmic censorship hypothesis 
and black-hole no-hair theorems, along with that of cuiTent 
open issues in relativistic astrophysics research, such as the 
mechanism responsible for gamma-ray bursts. Furthermore, 
numerical simulations of stellar gravitational collapse to black 
holes provide a unique mean of computing the gravitational 
waveforms emitted in such events, believed to be among the 
most important sources of detectable gravitational radiation. 

However, the modelling of black-hole spacetimes with col- 
lapsing matter-sources in multidimensions is one of the most 
formidable efforts of numerical relativity. This is due, on one 
hand, to the inherent difficulties and complexities of the sys- 
tem of equations which is to be integrated, the Einstein field 
equations coupled to the general-relativistic hydrodynamics 
equations, and, on the other hand, to the immense computa- 
tional resources needed to integrate the equations in the case 
of three-dimensional (3D) evolutions. In addition to the prac- 
tical difficulties encountered in the accurate treatment of the 
hydrodynamics involved in the gravitational collapse of a ro- 
tating neutron star to a black hole, the precise numerical com- 
putation of the gravitational radiation emitted in the process is 
particularly challenging as the energy released in gravitational 



waves is much smaller than the total rest-mass energy of the 
system. 

The ability to perform long-term numerical simulations of 
self-gravitating systems in general relativity strongly depends 
on the formulation adopted for the Einstein equations. The co- 
variant nature of these equations (the "many-fingered time" of 
relativity) leads to difficulties in constructing an appropriate 
coordinate representation which would allow for stable and 
accurate simulations. Over the years, the standard approach 
has been mainly based upon the unconstrained solution of the 
3+1 ADM formulation of the field equations, which, despite 
large-scale and dedicated collaborations 11,2,3] has gradually 
been shown to lack the stability properties necessary for long- 
term numerical simulations. In recent years, however, a con- 
siderable effort has been invested in extending the set of ADM 
equations solved by including at some level the solution of the 
constraint equations on each spatial hypersurface or 
by reformulating the ADM approach in order to achieve long- 
term stability (see, e.g., 0] and references therein). Building 
on the experience developed with lower-dimensional formula- 
tions, Nakamura, Oohara and Kojima ^ presented in 1987 a 
conformal traceless reformulation of the ADM s yste m which 
subsequent authors (see, e.g., 1 9 , To', T7, "12, T?, 14, 1?,'!?]) 
gradually showed to be robust enough to accomplish such a 
goal for different classes of spacetimes, including black holes 
and neutron stars (both isolated and in coalescing binary sys- 
tems). The most widespread version developed from this for- 
malism, which we refer to here as the NOK formulation, was 
given by fioll and is commonly referred to as the BSSN 
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formulation. 

In addition to the improvements achieved in the formulation 
of the field equations, successful long-term 3D evolutions of 
black holes in vacuum have been obtained in the last few years 
using excision techniques (see, e.g., lil7i..l8..19.,.2Qi.i2 uEzEsl 
I24ll2a.l23l '). although the original idea is much older l27ll . In 
this approach, the spacetime region within the black-hole hori- 
zon, and so causally disconnected, can be safely ignored with- 
out affecting the evolution outside the horizon as long as suit- 
able boundary conditions are specified at the excision surface. 
The simulations presented here show the applicability of exci- 
sion techniques also in non-vacuum spacetimes, namely dur- 
ing the collapse of rotating neutron stars to Kerr black holes. 
The excision technique, which is applied once the black-hole 
apparent horizon is found, permits to extend considerably the 
lifetime of the simulations, at least with the resolutions used 
here. This, in turn, allows for an accurate investigation of the 
dynamics of the trapped surfaces formed during the collapse, 
from which important information on the mass and spin of the 
black hole, as well as on the amount of energy which is lost to 
gravitational radiation, can be extracted. While our study was 
nearing completion, we have learnt that a similar approach has 
also been implemented with success f2§R . 

The presence of rotation in the collapsing stellar models 
requires multidimensional investigations, either in axisym- 
metry or in full 3D. The numerical investigations of black- 
hole formation (beyond spherical symmetry) started in the 
early 1980's with the pioneering work of Nakamura |29]. He 
adopted the (2h-1)h-1 formulation of the Einstein equations in 
cylindrical coordinates and introduced regularity conditions to 
avoid divergences at coordinate singularities. Nakamura used 
a "hypergeometric" slicing condition which prevents the grid 
points from reaching the singularity when a black hole forms. 
The simulations could track the evolution of the collapse of a 
IOMq "core" of a massive star with different amounts of ro- 
tational energy, up to the formation of a rotating black hole. 
However, the numerical scheme employed was not accurate 
enough to compute the emitted gravitational radiation. In 
subsequent works, Nakamura |30] (see also |8]) considered 
a configuration consisting of a neutron star with an accreting 
envelope, which was thought to mimic mass fall-back in a su- 
pernova explosion. Rotation and infall velocity were added 
to such a configuration, simulating an evolution dependent on 
the prescribed rotation rates and rotation laws. 

Later on, in a series of papers El mm HI, Bardeen, 

Stark and Piran studied the collapse of rotating relativistic 
polytropes to black holes, succeeding in computing the as- 
sociated gravitational radiation. The gravitational field and 
hydrodynamics equations were formulated using the 3 + 1 
formalism in two spatial dimensions, using the radial gauge 
and a mixture of singularity-avoiding polar and maximal slic- 
ings. The initial model was a spherically symmetric relativis- 
tic polytrope of mass AI in equilibrium. The gravitational 
collapse was induced by lowering the pressure in the initial 
model by a prescribed (and often very large) fraction and by 
simultaneously adding an angular momentum distribution ap- 
proximating rigid-body rotation. Their simulations showed 
that Kerr black-hole formation occurs only for angular mo- 



menta less than a critical value. Furthermore, the energy AE 
carried away through gravitational waves from the collapse to 
a Ken- black hole was found to be AE/Mc^ < 7 x 10"*, 
the shape of the waveforms being nearly independent of the 
details of the collapse. 

The axisymmetric codes employed in the aforementioned 
works adopted curvilinear coordinate systems which may lead 
to long-term numerical instabilities at coordinate singulari- 
ties. These coordinate problems have deterred researchers 
from building 3D codes in spherical coordinates. Recently, 
a general-purpose method (called "cartoon"), has been pro- 
posed to enforce axisymmetry in numerical codes based on 
Cartesian coordinates and which does not suffer from stability 
problems at coordinate singularities |35|. It should be noted, 
however, that the stability properties of the cartoon method 
are not fully understood yet, as discussed by |36]. Using 
this method, Shibata |37] investigated the effects of rotation 
on the criterion for prompt adiabatic collapse of rigidly and 
differentially rotating polytropes to a black hole, finding that 
the criterion for black-hole formation depends strongly on the 
amount of angular momentum, but only weakly on its (initial) 
distribution. The effects of shock heating when using a non- 
isentropic equation of state (EOS hereafter) are important in 
preventing prompt collapse to black holes in the case of large 
rotation rates. 

More recently, Shibata has performed axisymmet- 

ric simulations of the collapse of rotating supramassive neu- 
tron stars to black holes for a wide range of polytropic EOSs 
and with an improved implementation of the hydrodynamics 
solver (based on approximate Riemann solvers) with respect 
to the original implementation used in JstIi . Parameterizing 
the "stiffness" of the EOS through the polytropic index N, 
the final state of the collapse is a Kerr black hole without 
any noticeable disc formation, when the polytropic index N 
is in the range 2/3 < TV < 2. Based on the specific an- 
gular momentum distribution in the initial star, Shibata has 
estimated an upper limit to the mass of a possible disc as be- 
ing less than 10""^ of the initial stellar mass fWj. Unfortu- 
nately, such small masses cannot currently be confirmed with 
the presently-available resolutions in 3D simulations on uni- 
form grids. 

Three-dimensional, fully relativistic simulations of the col- 
lapse of supramassive uniformly rotating neutron stars to ro- 
tating black holes were presented in 1, 1 1.1 . The simulations fo- 
cused on = 1 polytropes and showed no evidence of mas- 
sive disc formation or outflows. These results are in agree- 
ment with those obtained in axisymmetry 1 12, 381 and with 
the new simulations reported by E3l (both in axisymmetry 
and in 3D) which show that for a rapidly rotating polytrope 
with J /AP < 0.9 (J being the angular momentum) all the 
mass falls promptly into the black hole, with no disc be- 
ing formed. Hence, all existing simulations agree that mas- 
sive disc formation from the collapse of neutron stars re- 
quires differential rotation, at least for a polytropic EOS with 
1 < TV < 2. 

Here, we present the results of new, fully 3D simulations 
of gravitational collapse of uniformly rotating neutron stars, 
both secularly and dynamically unstable, which we model as 
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relativistic polytropes. The angular velocities of our sam- 
ple of initial models range from slow rotation to the mass- 
shedding limit. For the first time in such 3D simulations, 
we have detected the event horizon of the forming black hole 
and showed that it can be used to achieve a more accurate 
determination of the black-hole mass and spin than it would 
be otherwise possible using the area of the apparent horizon. 
We have also considered several other approaches to measure 
the properties of the newly formed Kerr black hole, including 
the recently proposed isolated and dynamical-horizon frame- 
works. A comparison among the different methods has indi- 
cated that the dynamical-horizon approach is simple to imple- 
ment and yields estimates which are accurate and more robust 
than those of the equivalent methods. 

The simulations are performed with a new general relativis- 
tic hydrodynamics code, the Whisky code, in which the Ein- 
stein and hydrodynamics equations are finite-differenced on 
a Cartesian grid and solved using state-of-the-art numerical 
schemes (a first description of the code was given in |39]). 
The code incorporates the expertise developed over the past 
few years in the numerical solution of the Einstein equations 
and of the hydrodynamics equations in a curved spacetime 
(see ifijl Hill , but also |40] and references therein) and is 
the result of a collaboration among several European Insti- 
tutes |41]. 

As mentioned before, we have implemented in the Whisky 
code a robust excision algorithm which warrants the exten- 
sion of the lifetime of the simulations far beyond the evolu- 
tion times when the black holes first form. Our calculations 
are starting from initially axisymmetric stellar models but are 
performed in full 3D to allow for departures from the initial 
axial symmetry. Overall, our results show that the dynamics 
of the collapsing matter is strongly influenced by the initial 
amount of angular momentum in the progenitor neutron star, 
which, when sufficiently high, leads to the formation of an 
unstable flattened object. All of the simulations performed for 
realistic initial data and a polytropic equation of state show no 
evidence of shock formation preventing a prompt collapse to a 
black hole, nor the presence of matter on stable orbits outside 
the black hole. It should be remarked, however, that both of 
these conclusions may change if the initial stellar models are 
rotating differentially. 

The use of numerical grids with uniform spacing and 
the present computational resources have placed the outer 
boundary of our computational box in regions of the space- 
time where the gravitational waves have not yet reached 
their asymptotic form. As a result, the information on the 
gravitational waveforms that we extract through perturbative 
techniques does not provide interesting information 

besides the obvious change in the stellar quadrupole mo- 
ment. Work is now in progress to use mesh refinement tech- 
niques 1 44] to move the outer boundary sufficiently far from 
the source so that important information can be extracted on 
the gravitational wave emission produced during the collapse. 
The results of these investigations will be presented in a com- 
panion paper i4^ . 

The paper is organized as follows: Section HI] describes 
the formulation we adopt for the Einstein and hydrodynam- 



ics equations, the way they are implemented in the code 
and a brief discussion of how the excision techniques can be 
employed within a hydrodynamical treatment making use of 
high-resolution shock-capturing (HRSC) schemes. To avoid 
detracting the reader's attention from the physical problem 
considered here, we have confined most of the technical de- 
tails concerning the numerical implementation of the hydro- 
dynamical equations to Appendix |X| Section [111] is therefore 
devoted to describing the various properties of the initial stel- 
lar models. The following two Sections, llVl andlVl present 
our results regarding the dynamics of the collapsing stars and 
of the trapped surfaces, respectively. In both cases we will 
be considering and comparing the dynamics of slowly and 
of rapidly rotating stellar models. The paper ends with Sec- 
tion |V1J which contains a summary of the results obtained 
and the perspectives for further investigations. Finally, Ap- 
pendix |B] is devoted to presenting some of the tests the code 
passes with very high accuracy. 

We here use a spacelike signature (— , +, +, +) and a sys- 
tem of units in which c = G = Mq = 1 (unless explicitly 
shown otherwise for convenience). Greek indices are taken to 
run from to 3, Latin indices from 1 to 3 and we adopt the 
standard convention for the summation over repeated indices. 



II. BASIC EQUATIONS AND THEIR IMPLEMENTATION 

The Whisky code solves the general relativistic hydrody- 
namics equations on a 3D numerical grid with Cartesian co- 
ordinates. The code has been constructed within the frame- 
work of the Cactus Computational Toolkit (see |46] for de- 
tails), developed at the Albert Einstein Institute (Golm) and 
at the Louisiana State University (Baton Rouge). This pub- 
lic domain code provides high-level facilities such as paral- 
lelization, input/output, portability on different platforms and 
several evolution schemes to solve general systems of partial 
differential equations. Clearly, special attention is dedicated 
to the solution of the Einstein equations, whose matter-terms 
in non-vacuum spacetimes are handled by the Whisky code. 
While the Whisky code is entirely new, its initial develop- 
ment has benefitted in part from the release of a public ver- 
sion of the general relativistic hydrodynamics code described 
in Il4[|47ll . and developed mostly by the group at the Wash- 
ington University (St. Louis). 

The Whisky code, however, incorporates important re- 
cent developments regarding, in particular, new numerical 
methods for the solution of the hydrodynamics equations that 
have been described in detail in |39| and will be briefly 
reviewed in Appendix |X| These include: ( i) the Piece- 
wise ParaboHc Method (PPM ) |.4H and the Essentially Non- 
Oscillatory (ENO) methods 14^ for the cell-reconstruction 
procedure; ( ii) the Harten-Lax-van Leer-Einfeldt (HLLE) 1 53] 
approximate Riemann solver, the Marquina flux formiila ] 5 ij ; 
( Hi) the analytic expression for the left eigenvectors and 
the compact flux formulae 1.5 3il for a Roe-type Riemann solver 
and a Marquina flux formula; (iv) the use of a "method of 
lines" (MoL) approach for the implementation of high-order 
time evolution schemes; (v) the possibility to couple the gen- 
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eral relativistic hydrodynamics equations with a conformally 
decomposed three-metric. The incorporation of these new nu- 
merical techniques in the code has led to a much improved 
ability to simulate relativistic stars, as shown in Appendix 151 
which is devoted to code tests. 

While the Cactus code provides at each time step a solu- 
tion of the Einstein equations flsil 



(2.1) 



where G^i, is the Einstein tensor and T^^ is the stress-energy 
tensor, the Whisky code provides the time evolution of the 
hydrodynamics equations, expressed through the conservation 
equations for the stress-energy tensor T^'' and for the matter 
current density 



0, 



0. 



(2.2) 



In what follows we briefly discuss how both the right and 
the left-hand side of equations il.H are computed within the 
coupled Cactus/Whisky codes. 



A. Evolution of the field equations 

We here give only a brief overview of the system of equa- 
tions for the evolution of the field equations, but refer the 
reader to for more details. Many different formulations of 
the equations have been proposed throughout the years, start- 
ing with the ADM formulation in 1962 js^l- As mentioned 
in the Introduction, we use the NOK f§\ formulation, which 
is based on the ADM construction and has been further devel- 
oped in |9]. 

In the ADM formulation fs?], the spacetime is foliated with 
a set of non-intersecting spacelike hypersurfaces. Two kine- 
matic variables relate the hypersurfaces: the lapse function a, 
which describes the rate of advance of time along a timelike 
unit vector normal to a spacelike hypersurface, and the 
shift three-vector /3* that relates the coordinates of two space- 
like hypersurfaces. In this construction the fine element reads 



-{a^ - I3il3')dt^ + 2l3,dx'dt + jijdx'dx^ . (2.3) 



The original ADM formulation casts the Einstein equations 
into a first-order (in time) quasi-Unear system of equa- 
tions. The dependent variables are the three-metric 7^ and 
the extrinsic curvature Kij, with first-order evolution equa- 
tions given by 



^a^J = -2aK,j +V,(3j +Vj(3, 



(2.4) 



dtK, 



-Stt S, 



Ri 



K K, . 



IK K"' 



Here, denotes the covariant derivative with respect to the 
three-metric 7^ , i?^ is the Ricci curvature of the three-metric, 
K = 'j'^^ Kij is the trace of the extrinsic curvature, Sij is the 
projection of the stress-energy tensor onto the spacelike hy- 
persurfaces and S = 7'-' Sij (for a more detailed discussion, 
see isfiil '). In addition to the evolution equations, the Einstein 
equations also provide four constraint equations to be satisfied 
on each spacelike hypersurface. These are the Hamiltonian 
constraint equation 



167rp^ 



, 



and the momentum constraint equations 



Snf = . 



(2.6) 



(2.7) 



(2.5) 



In equations (l2.4> -( l277t . p^dm ™d f ^re the energy density 
and the momentum density as measured by an observer mov- 
ing orthogonally to the spacelike hypersurfaces. 

Details of our particular implementation of the confor- 
mal traceless reformulation of the ADM system as proposed 
by L8.. „9„ .1Q.1 are extensively described in 1 13. ,57.1 and will not 
be repeated here. We only mention, however, that this formu- 
lation makes use of a conformal decomposition of the three- 
metric, jij = e~^'^7ij, and the trace-free part of the extrinsic 
curvature, Aij — Kij — jijK/3, with the conformal factor (p 
chosen to satisfy e'*'^ = 7^^'^, where 7 is the determinant of 
the spatial three-metric 7^ . In this formulation, in addition to 
the evolution equations for the conformal three-metric 7^ and 
the conformal traceless extrinsic curvature Aij, there are evo- 
lution equations for the conformal factor (f>, for the trace of the 
extrinsic curvature K and for the "conformal connection func- 
tions" = 7*-' j . We note that although the final mixed, first 

and second-order, evolution system for if, 7^ , Jiy , f *| 

is not in any immediate sense hyperbolic, there is evidence 
showing that the formulation is at least equivalent to a hy- 
perbolic system l,5&. i59t i6Qi] . In the formulation of |9], the 
auxiliary variables Fi = — 7^ j were used instead of the 

V\ 

In Refs. fT3','6l'l the improved properties of this conformal 
traceless formulation of the Einstein equations were compared 
to the ADM system. In particular, in 1 1 3 ] a number of strongly 
gravitating systems were analysed numerically with conver- 
gent HRSC methods with total-variation-diminishing (TVD) 
schemes using the equations described in |47]. These included 
weak and strong gravitational waves, black holes, boson stars 
and relativistic stars. The results showed that this treatment 
leads to numerical evolutions of the various strongly gravitat- 
ing systems which did not show signs of numerical instabil- 
ities for sufficiently long times. However, it was also found 
that the conformal traceless formulation requires grid resolu- 
tions higher than the ones needed in the ADM formulation to 
achieve the same accuracy, when the foliation is made using 
the "if-driver" approach discussed in |62]. Because in long- 
term evolutions a small error-growth rate is the most desirable 
property, we have adopted the conformal traceless formulation 
as our standard form for the evolution of the field equations. 
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1. Gauge choices 

The code is designed to handle arbitrary shift and lapse con- 
ditions, which can be chosen as appropriate for a given space- 
time simulation. More information about the possible families 
of spacetime slicings which have been tested and used with 
the present code can be found in 1 13l l22ll . Here, we limit our- 
selves to recalling details about the specific foUations used in 
the present evolutions. In particular, we have used hyperbolic 
if-driver slicing conditions of the form 

{dt-l3'di)a^~f{a)a^{K~Ko), (2.8) 

with f{a) > and = K{t = 0). This is a generalization 
of many well known slicing conditions. For example, setting 
/ = 1 we recover the "harmonic" sUcing condition, while, 
by setting / = q/a, with q an integer, we recover the gener- 
alized "1+log" slicing condition |^2J- In particular, all of the 
simulations discussed in this paper are done using condition 
(12. 8> with / = 2/ a. This choice has been made mostly be- 
cause of its computational efficiency, but we are aware that 
"gau ge p athologies" could develop with the "1+log" slic- 
ings |64("65'l. 

As for the spatial gauge, we use one of the "Gamma-driver" 
shift conditions proposed in |22] (see also |57|), that essen- 
tially act so as to drive the to be constant. In this re- 
spect, the "Gamma-driver" shift conditions are similar to the 
"Gamma-freezing" condition Stf^ — 0, which, in turn, is 
closely related to the well-known minimal distortion shift con- 
dition ll66ll . The differences between these two conditions in- 
volve the Christoffel symbols and are basically due to the fact 
that the minimal distortion condition is covariant, while the 
Gamma-freezing condition is not. 

In particular, all of the results reported here have been ob- 
tained using the hyperbolic Gamma-driver condition, 

dfP' ^Fdtt' -r^dtP\ (2.9) 

where F and -q are, in general, positive functions of space 
and time. For the hyperbolic Gamma-driver conditions it is 
crucial to add a dissipation term with coefficient -q to avoid 
strong oscillations in the shift. Experience has shown that by 
tuning the value of this dissipation coefficient it is possible to 
almost freeze the evolution of the system at late times. We 
typically choose F = | and rj — 'i and do not vary them in 
time. 



B. Evolution of the hydrodynamics equations 

An important feature of the Whisky code is the imple- 
mentation of a conservative formulation of the hydrodynamics 
equations |52, 67, 68|, in which the set of equations (12. 2> is 
written in a hyperbolic, first-order and flux-conservative form 
of the type 

9tq + 5.fW(q) -s(q) , (2.10) 

where f (q) and s(q) are the flux-vectors and source terms, 
respectively L40>l . Note that the right-hand side (the source 



terms) depends only on the metric, and its first derivatives, 
and on the stress-energy tensor. Furthermore, while the sys- 
tem ( I2.10> is not strictly hyperbolic, strong hyperbolicity is 
recovered in a flat spacetime, where s(q) = 0. 

As shown by 1 68 ] , in order to write system ( 12. 2> in the form 
of system (I2.10> . the primitive hydrodynamical variables {i.e. 
the rest-mass density p and the pressure p (measured in the 
rest-frame of the fluid), the fluid three-velocity (measured 
by a local zero-angular momentum observer), the specific in- 
ternal energy e and the Lorentz factor W) are mapped to the 
so called conserved variables q = {D,S'',t) via the relations 

D = ^Wp, 

S' = ^phW^v' , (2.11) 
T = ^{phW^~p)^D, 

where h = l + e + p/pis the specific enthalpy and 
W = {\ ~ ^ijV^v^)^^/^ . Note that only five of the seven 
primitive variables are independent. 

In order to close the system of equations for the hydrody- 
namics an EOS which relates the pressure to the rest-mass 
density and to the energy density must be specified. The code 
has been written to use any EOS, but all of the tests so far have 
been performed using either an (isentropic) polytropic EOS 

p ^ Kp^ , (2.12) 

e = P+Y&i^ (2-13) 

or an "ideal fluid" EOS 

p=(r-l)pe. (2.14) 

Here, e is the energy density in the rest-frame of the fluid, 
K the polytropic constant (not to be confused with the trace 
of the extrinsic curvature defined earlier) and F the adia- 
batic exponent. In the case of the polytropic EOS ( I2.12t . 
r = 1 + where N is the polytropic index and the 

evolution equation for r needs not be solved. In the case of 
the ideal-fluid EOS ( 12.141 1. on the other hand, non-isentropic 
changes can take place in the fluid and the evolution equa- 
tion for T needs to be solved. In addition to the EOSs ( I2.12t 
and ( 12.141 1. a "hybrid" EOS (suitable for core-collapse simu- 
lations), as described in 1 69, 70], has also been implemented. 

Note that polytropic EOSs ( I2.12t do not allow any transfer 
of kinetic energy to thermal energy, a process which occurs in 
physical shocks (shock heating). However, we have verified, 
by performing simulations with the more general EOS ( I2.14t 
on some selected cases, that for the physical systems treated 
here, shock heating is not important (no shocks form during 
our simulations). Since in our numerical scheme using e gen- 
eral EOS like ( I2.14t is more expensive than using a polytropic 
EOS, the systematic investigations presented here have been 
obtained using the latter. 

Additional details of the formulation we use for the hydro- 
dynamics equations can be found in i40ll . We stress that an 
important feature of this formulation is that it has allowed to 
extend to a general relativistic context the powerful numer- 
ical methods developed in classical hydrodynamics, in par- 
ticular HRSC schemes based on linearized Riemann solvers 
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(see Such schemes are essential for a correct represen- 

tation of shocks, whose presence is expected in several as- 
trophysical scenarios. Two important results corroborate this 
view. The first one, by Lax and Wendroff |71], states that 
a stable scheme converges to a weak solution of the hydrody- 
namical equations. The second one, by Hou and LeFloch 1731 . 
states that, in general, a non-conservative scheme will con- 
verge to the wrong weak solution in the presence of a shock, 
hence underlining the importance of flux-conservative formu- 
lations. For a full introduction to HRSC methods the reader is 
also referred to 



C. Hydrodynamical excision 

Excision boundaries are usually based on the principle that 
a region of spacetime that is causally disconnected can be ig- 
nored without this affecting the solution in the remaining por- 
tion of the spacetime. Although this is true for signals and 
perturbations travelling at physical speeds, numerical calcu- 
lations may violate this assumption and disturbances, such 
as gauge waves, may travel at larger speeds thus leaving the 
physically disconnected regions. 

Note that, in non-vacuum spacetimes, the excision bound- 
aries for the hydrodynamical and the metric fields need not 
be the same. For the fluid quantities, in fact, all character- 
istics emanating from an event in spacetime will propagate 
within the sound-cone at that event, and, for physically real- 
istic EOSs, this sound-cone will always be contained within 
the light-cone at that event. As a result, if a region of space- 
time contains trapped surfaces, both the hydrodynamical and 
the metric fields are causally disconnected and both can be ex- 
cised there. On the other hand, there may be situations, e.g., 
when the bulk flow is locally supersonic but no trapped sur- 
faces are present, in which it is possible (at least in principle) 
to excise the hydrodynamical fields without having to do the 
same for the metric fields. We have not used this option here 
and the hydrodynamical excision implemented in our simula- 
tions has always been made within regions of the spacetime 
contained inside trapped surfaces. 

A first naive implementation of an excision algorithm 
within a HRSC method could ensure that the data used to con- 
struct the flux at the excision boundary is extrapolated from 
data outside the excision region. This may appear to be a 
good idea since HRSC methods naturally change the stencils 
depending on the data locally. In general, however, this ap- 
proach is not guaranteed to reduce the total variation of the 
solution and simple examples may be produced that fail with 
this boundary condition. 

An effective solution, however, is not much more compli- 
cated and can be obtained by applying at the excision bound- 
ary the simplest outflow boundary condition (here, by outflow 
we mean flow into the excision region). In practice, we apply 
a zeroth-order extrapolation to all variables at the boundary, 
i.e. a simple copy of the hydrodynamical variables across the 
excision boundary (note that setting the hydrodynamical fields 
inside the excised region to zero would still yield an outflow 
boundary condition, but leads to incorrect outflow speeds). If 



the reconstruction method requires more cells inside the ex- 
cision region, we force the stencil to only consider the data 
in the exterior and the first interior cell. Although the actual 
implementation of this excision technique may depend on the 
reconstruction method used, the working principle is always 
the same. 

The location of the excision boundary itself is based on 
the determination of the apparent horizon which, within the 
Cactus code, is obtained using the fast finder of Thornburg 
|l76l|. The excision boundary is placed a few gridpoints (typi- 
cally 4), within a surface which is 0.6 times the size of the ap- 
parent horizon. This may not be a suitable outflow boundary 
on a Cartesian grid, as pointed out by L77.,78.1 . However, sim- 
ilar or larger excision regions show no problems in vacuum 
evolutions and since the sound-cones are always contained 
within the light-cones, we expect no additional problems to 
arise from the hydrodynamics. 

More details on how the hydrodynamical excision is ap- 
plied in practice, as well as tests showing that this method is 
stable, consistent and converges to the expected order will be 
published in a separate paper . 



III. INITIAL STELLAR MODELS 

As mentioned earlier, this paper is specially dedicated to 
study the gravitational collapse of slowly and rapidly rotat- 
ing supramassive relativistic stars, in uniform rotation, that 
have become unstable to axisymmetric perturbations. Given 
equilibrium models of gravitational mass M and central en- 
ergy density Cc along a sequence of fixed angular momen- 
tum or fixed rest mass, the Friedman, Ipser & Sorkin crite- 
rion dM/dec — 1 80] can be used to locate the exact on- 
set of the secular instability to axisymmetric collapse. The 
onset of the dynamical instability to collapse is located near 
that of the secular instability but at somewhat larger central 
energy densities. Unfortunately, no simple criterion exists to 
determine this location, but the expectation mentioned above 
has been confirmed by the simulations performed here and by 
those discussed in lllll . Note that in the absence of viscosity 
or strong magnetic fields, the star is not constrained to rotate 
uniformly after the onset of the secular instability and could 
develop differential rotation. In a realistic neutron star, how- 
ever, viscosity or intense magnetic fields are likely to enforce 
a uniform rotation and cause the star to collapse soon after it 
passes the secular instability limit. 

The initial data for our simulations are constructed using a 
2D numerical code, that computes accurate stationary equi- 
librium solutions for axisymmetric and rapidly rotating rel- 
ativistic stars in polar coordinates |81]. The data are then 
transformed to Cartesian coordinates using standard coordi- 
nate transformations. The same initial data routines have been 
used in previous 3D simulations 1 13, 14, 82] and details on the 
accuracy of the code can be found in ] 83]. 

For simplicity, we have focused on initial models con- 
structed with the polytropic EOS ( I2.12> . choosing F = 2 and 
polytropic constant Kid — 100 to produce stellar models that 
are, at least qualitatively, representative of what is expected 
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FIG. 1: Gravitational mass shown as a function of the central energy 
density for equilibrium models constructed with the polytropic EOS, 
for r = 2 and polytropic constant Kjo = 100. The solid, dashed 
and dotted lines correspond to the sequence of non-rotating models, 
the sequence of models rotating at the mass-shedding limit and the 
sequence of models that are at the onset of the secular instability 
to axisymmetric perturbations. Also shown are the secularly (open 
circles) and dynamically unstable (filled circles) initial models used 
in the collapse simulations. 



from observations of neutron stars. More specifically, we have 
selected four models located on the line defining the onset of 
the secular instability and having polar-to-equatorial axes ra- 
tio of roughly 0.95, 0.85, 0.75 and 0.65 (these models are indi- 
cated as S1-S4 in FigQ, respectively. Four additional models 
were defined by increasing the central energy density of the 
secularly unstable models by 5%, keeping the same axis ratio. 
These models (indicated as D1-D4 in FigQ were expected 
and have been found to be dynamically unstable. 

Figure [2 shows the gravitational mass as a function of 
the central energy density for equilibrium models constructed 
with the chosen polytropic EOS. The solid, dashed and dot- 
ted lines correspond respectively to: the sequence of non- 
rotating models, the sequence of models rotating at the mass- 
shedding limit and the sequence of models that are at the onset 
of the secular instability to axisymmetric perturbations. Fur- 
thermore, the secularly and dynamically unstable initial mod- 
els used in the collapse simulations are shown as open and 
filled circles, respectively. 

Table m summarizes the main equilibrium properties of the 
initial models. The circumferential equatorial radius is de- 
noted as i?e, while 51 is the angular velocity with respect 
to an inertial observer at infinity, and rp/r^ is the ratio of 
the polar to equatorial coordinate radii. The height of the 
corotating innermost stable circular orbit (ISCO) is defined 
as /i+ = i?+ — i?e, where _R+ is the circumferential radius 



TABLE I: Equilibrium properties of the initial stellar models. The 
different columns refer respectively to: the central rest-mass density 
Pc the ratio of the polar to equatorial coordinate radii rp/r^, the 
gravitational mass A/, the circumferential equatorial radius Re, the 
angular velocity Q, the ratio J/M^ where J is the angular momen- 
tum, the ratio of rotational kinetic energy to gravitational binding 
energy T/|VF|, and the "height" of the corotating ISCO h+. All 
models have been computed with a polytropic EOS with A'id = 100 
and r = 2. 



Model 


Pl 


Vp/re 


M 


Re 


j/M^ r/|w|t 


h+ 


SI 


3.154 


0.95 


1.666 7.82 1.69 


0.207 


1.16 


1.18 


S2 


3.066 


0.85 


1.729 


8.30 2.83 


0.363 


3.53 


0.51 


S3 


3.013 


0.75 


1.798 


8.90 3.49 


0.470 


5.82 


0.04 


S4 


2.995 


0.65 


1.863 


9.76 3.88 


0.545 


7.72 




Dl 


3.280 


0.95 


1.665 


7.74 1.73 


0.206 


1.16 


1.26 


D2 


3.189 


0.85 


1.728 


8.21 2.88 


0.362 


3.52 


0.58 


D3 


3.134 


0.75 


1.797 


8.80 3.55 


0.468 


5.79 


0.10 


D4 


3.116 


0.65 


1.861 


9.65 3.95 


0.543 


7.67 





t X lO-^* 
t X 10"^ 



for a corotating ISCO observer Note that in those models for 
which a value of is not reported, all circular geodesic orbits 
outside the stellar surface are stable. Other quantities shown 
are the central rest-mass density pc, the angular momentum 
J, and the ratio of rotational kinetic energy to gravitational 
binding energy r/|M^|. 



IV. DYNAMICS OF THE MATTER 

This Section discusses the dynamics of the matter during 
the collapse of the initial stellar models described in the pre- 
ceding section. All of the simulations reported here have been 
computed using a uniformly spaced computational grid for 
which symmetry conditions are imposed across the equatorial 
plane. Different spatial resolutions have been used to check 
convergence and improve the accuracy of the results, with the 
finest resolution having been obtained using 288'^ x 144 cells. 
While the precise numbers depend on the resolution used and 
on the model simulated, as a general rule we have used ^ 50% 
of the cells in the x-direction to cover the star in case Dl and 
~ 66% of the cells in the x-direction to cover the star in D4. 
As a result, the outer boundary is set at ~ 2.0 times the stellar 
equatorial radius for Dl and at ^ 1.4 times the stellar equato- 
rial radius for D4. 

The hydrodynamics equations have been solved employing 
the Marquina flux formula and a third-order PPM reconstruc- 
tion, which was shown in |84] to be superior to other meth- 
ods in maintai ning a highly-accurate angular-velocity profile 
(see also 1 8^ l85ll for recent 3D evolutions of rotating rel- 
ativistic stars with a third-order order PPM reconstruction). 
The Einstein field equations, on the other hand, have been 
evolved using the ICN evolution scheme, the "1+log" slic- 
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FIG. 2: 1/2 norm of the Hamiltonian constraint violation for the ini- 
tial model Dl shown as a function of time. The different lines refer 
to different grid resolutions, but in all cases the IVP was solved after 
the pressure was uniformly decreased to trigger the collapse. 



ing condition and the "Gamma-driver" shift conditions fl^. 
Finally, both polytropic and ideal-fluid EOSs have been used, 
although no significant difference has been found in the dy- 
namics between the two cases. This is most probably related 
to the small J/AP of the uniformly rotating initial models 
considered here. This implies a relatively rapid collapse and 
as a result we do not see any shock formation (see below for 
a more complete discussion). Hereafter we will restrict our 
attention to a polytropic EOS only. 

Given an initial stellar model which is dynamically unsta- 
ble, simple round-off errors would be sufficient to produce an 
evolution leading either to the gravitational collapse to a black 
hole or to the migration to the stable branch of the equilibrium 
configurations I14ll (we recall that both evolutions are equally 
probable mathematically, although it is easier to imagine that 
it would collapse to a black hole). In general, however, leaving 
the onset of the dynamical instability to the cumulative effect 
of the numerical truncation eiTor is not a good idea, since this 
produces instability-growth times that are dependent on the 
grid-resolution used. 

For this reason, we induce the collapse by slightly reduc- 
ing the pressure in the initial configuration. This is done uni- 
formly throughout the star by using a polytropic constant for 
the evolution K that is smaller than the one used to calculate 
the initial data A'id- The accuracy of the code is such that only 
very small perturbations are sufficient to produce the collapse 
and we have usually adopted (i^iD — K)/Kn:, < 2%. 

After imposing the pressure reduction, the Hamiltonian and 
momentum constraints are solved to enforce that the con- 
straint violation is at the truncation-error level. We refer to 



this procedure as to solving the initial value problem (IVP), 
which ensures that second order convergence holds from the 
start of the simulations, as shown in Fig. |2] for the L2 norm 
of the Hamiltonian constraint. Strict second-order conver- 
gence is lost when the excision is introduced, although the 
code remains convergent at a lower rate while the norms of 
the Hamiltonian constraint start to grow exponentially (this is 
not shown in Fig.|2j- We are presently investigating the origin 
of the deterioration of the convergence rate at the time of ex- 
cision, although this is somewhat unavoidable when excising 
a spherical region in a Cartesian rectangular grid in the course 
of the evolution. 

Details on how we solve the IVP implementing the York- 
Lichnerowicz conformal transverse-traceless decomposition 
can be found in [86]. If, on the other hand, the IVP is not 
solved after the pressure change, the constraints violation in- 
creases twice as fast and converges to second order only after 
an initial period of about 20 M ^ 0.17 ms. To assess the 
validity of our procedure to trigger the collapse, we also per- 
form the pressure change after the evolution has started and 
without solving the IVP. In this case, after the system has re- 
covered from the perturbation, the violation of the constraints 
is only a few percent different from the case in which the IVP 
is solved. Furthermore, other dynamical features of the col- 
lapse, such as the instant at which the apparent horizon is first 
formed (see Section|V]for a detailed discussion), do not vary 
by more than 1%. 

The dynamics resulting from the collapse of models Sl- 
S4 and D1-D4 are extremely similar and no qualitative dif- 
ferences have been detected. However, as one would expect, 
models D1-D4 collapse more rapidly to a black hole (the for- 
mation of the apparent horizon appears about 5% earlier in co- 
ordinate time), are computationally less expensive and there- 
fore better suited for a detailed investigation. As a result, in 
what follows we will restrict our discussion to the collapse of 
the dynamically unstable models and distinguish the dynam- 
ics of case Dl, in Section HV Al from that of model D4, in 
Section llV^ 



A. Slowly rotating stellar models 

We start by discussing the dynamics of the matter by look- 
ing at the evolution of the initial stellar model Dl which is 
slowly rotating (thus almost spherical, with rp/r^ = 0.95) 
and has the largest central density (cf. Fig.[2and Table|l|l. 

We show in Figs.|3j0some representative snapshots of the 
evolution of this initial model. The different panels of Fig. 13 
in particular, refer to the initial and intermediate stages of the 
collapse and show the isocontours of the rest-mass density and 
velocity field in the (x, y) plane (left column) and in the {x, z) 
plane (right column), respectively. The isobaric surfaces are 
logarithmically spaced starting from p = 2.0 x 10^^ and go- 
ing up to p = 2.0 X 10^"^ at the stellar interior. The velocity 
vector field is expressed in units of c and the length for a ve- 
locity of 0.2 c is shown in the lower-right panel. Panels on the 
same row refer to the same instant in time and this is indicated 
in ms in the top-right corner of each panel. The units on both 
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FIG. 3: Collapse sequence for the slowly rotating model Dl. Different panels refer to different snapshots during the collapse and show the 
isocontours of the rest-mass density and velocity field in the {x, y) plane (left column) and in the [x, z) plane (right column), respectively. 
The isobaric surfaces are logarithmically spaced and a reference length for the vector field is shown in the lower right panel for a velocity of 
0.2 c. The time of the different snapshots appears in the top right corner of each panel and it is given in ms, while the units on both axes are 
expressed in km. See main text for a discussion and Fig.|5|for a comparison with the collapse of a rapidly rotating model. 
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FIG. 4: Magnified view of the final stages of tlie collapse of model Dl. A region around the singularity that has formed is excised from 
the computational domain and is shown as an area filled with squares. Also shown with a thick dashed line is the coordinate location of the 
apparent horizon. Note that because of rotation this surface is not a two-sphere, although the departures are not significant and cannot be 
appreciated from the Figure {cf. Fig. ll6l for a clearer view). 



axes are expressed in km and refer to coordinate lengths. This 
sequence has been obtained with a grid of 288^ x 144 zones 
but the data for the velocity field has been down-sampled to 
produce clearer figures. The data have also been restricted to 
a single octant in the (x, z) plane to provide a magnified view. 
In all cases, however, the whole extent of the numerical grid 
is reported in the figures. 

Overall, the sequence shown in Fig.|3]is simple to illustrate. 
During the collapse of model Dl spherical symmetry is almost 
preserved; as the star increases its compactness and the matter 
is compressed to larger pressures, the velocity field acquires a 
radial component (which was zero initially) that grows to rel- 
ativistic values. This is clearly shown in the panels at t = 0.49 
ms and t = 0.54 ms, which indicate that the star roughly pre- 
serves the ratio of its polar and equatorial radii (see also 
Fig.0, while radial velocities in excess of ^ 0.28 c can be 
easily reached. The behaviour of the angular velocity during 
this collapse will be analysed in more detail in Section ITV CI 
but we can here anticipate that it does not show appreciable 
departures from a profile which is uniform inside the star 

Soon after t = 0.54 ms, (i.e. at t = 0.546 ms = 66.6 M 
in the high-resolution run), an apparent horizon is found and 
when this has grown to a sufficiently large area, the portion 
of the computational domain containing the singularity is ex- 
cised. A discussion of how the trapped surfaces are studied 
in practice will be presented in Section |y] while details on 
the hydrodynamical excision have been given in Section HTd 
Here, we just remark that the use of an excised region and the 
removal of the singularity from the computational domain is 
essential for extending the calculations significantly past this 
point in time. Figure 0] shows a magnified view of the final 
stages of the collapse of model Dl. Indicated as an area filled 
with squares is the excised region of the computational do- 



main, which is an approximation of a sphere on the uniform 
Cartesian grid, i.e. a "lego-sphere". Also shown with a thick 
dashed line is the coordinate location of the apparent horizon 
and it should be remarked that, because of rotation, this sur- 
face is not a coordinate two-sphere, although the departures 
are not significant and cannot be appreciated in Fig. |4] (see 
Section|Vland TablelHlfor details), hit = 0.57 ms, the time 
which Fig.|4]refers to, most of the matter has already fallen 
within the apparent horizon and has assumed an oblate shape. 

The numerical calculations were carried out up to 
t ~ 0.73 ms ^ 89 M, thus using an excised region in a dy- 
namical spacetime for more than 26% of the total computing 
time. By this point, all of the stellar matter has collapsed well 
within the event horizon and the Hamiltonian constraint vio- 
lation has become very large. 

Overall, confirming what was already discussed by several 
authors in the past, the gravitational collapse of the slowly ro- 
tating stellar model Dl takes place in an almost spherical man- 
ner and we have found no evidence of shock formation which 
could prevent the prompt collapse to a black hole, nor appre- 
ciable deviations from axisymmetry {cf. left panel of Fig.|4}. It 
is possible, although not likely, that these qualitative features 
may be altered when a realistic EOS is used, since in this case 
shocks may appear, whose heating could stall or prevent the 
prompt collapse to a black hole. However, as mentioned in the 
Introduction, more dramatic changes are expected to appear if 
the initial configurations are chosen to have larger initial an- 
gular momenta and in particular when J/M^ > 1 fl2ll2^ . A 
first anticipation of the important corrections that centrifugal 
effects could produce is presented in the following Section, 
where we examine the dynamics of a rapidly rotating stellar 
model. 
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FIG. 5: Collapse sequence for the rapidly rotating model D4. The conventions used in these panels are the same as in Fig.|3| which can be 
used for a comparison with the collapse of a slowly rotating model. Note that a region around the singularity that has formed is excised from 
the computational domain and is indicated as an area filled with squares. Also shown with a thick dashed line is the coordinate location of the 
apparent horizon. 



B. Rapidly rotating stellar models 

We next consider the dynamics of the matter during the col- 
lapse of model D4 which, being rapidly rotating, is already 



rather flattened initially (i.e. rp/re — 0.65) and has the largest 
J/M'^ among the dynamically unstable models (cf. Fig.^and 
TableHJ. 

As for the slowly rotating star Dl, we show in Figs.|5l-|6l 
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FIG. 6: Magnified view of tiie final stages of tiie collapse of model D4. Note that a region around the singularity that has formed is excised 
from the computational domain and this is indicated as an area filled with squares. Also shown with a thick dashed line is the coordinate 
location of the apparent horizon. Note that because of the rapid rotation, this surface has significant departures from a two-sphere {cf. Fig. 1161 
for a clearer view). 



some representative snapshots of the evolution of this rapidly 
rotating model. The data has been computed using the same 
resolution of 288^ x 144 zones and the isocontour levels 
shown for the rest-mass density are the same used in Fig.|3}H 
It is apparent from the panels of Fig. [preferring to < = 0, that 
model D4 is considerably more oblate than Dl, as one would 
expect for a star rotating at almost the mass-shedding limit. 

Also in this case, we believe that the sequence in Fig.|5lis 
simple to illustrate. Note, in particular, how the dynamics is 
very similar to the one discussed for model Dl up to a time 
t ~ 0.49 ms. However, as the collapse proceeds, significant 
differences between the two models start to emerge and in the 
case of model D4 the large angular velocity of the progenitor 
stellar model produces significant deviations from a spheri- 
cal infall. Indeed, the parts of the star around the rotation 
axis that are experiencing smaller centrifugal forces collapse 
more promptly and, as a result, the configuration increases its 
oblateness. 

This is illustrated in Fig. which shows the time evolu- 
tion of the ratio of the polar and equatorial proper radii for all 
models in Table |I] (note that these ratios should not be con- 
fused with those in Table |l] that refer, instead, to coordinate 
radii). Each curve in Fig. extends until all of the matter 
along the z-axis has fallen inside the apparent horizon of the 
newly formed black hole. Clearly, in all cases the oblateness 
increases as the collapse proceeds and this is much more evi- 
dent for those stellar models that are rapidly rotating initially. 
In particular, for the most rapidly rotating models D4 and S4, 
the ratio between polar and equatorial proper radii becomes as 
small as 0.45 at the time when the matter on the rotation axis 
is below the apparent horizon. 

At about t = 0.64 ms (i.e. at t = 0.649 ms = 70.8 M in 
the high-resolution run), the collapse of model D4 produces 
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FIG. 7: Ratio of the proper polar radius to the proper equatorial ra- 
dius for all the initial models. Each curve ends at the time when, for 
each simulation, all the matter along the z-axis has fallen below the 
apparent horizon. 



an apparent horizon. Soon after this, the central regions of the 
computational domain are excised, preventing the code from 
crashing and thus allowing for an extended time evolution. 
The dynamics of the matter at this stage is shown in the lower 
panels of Fig.|5j which refer io t ~ 0.67 ms and where both 
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the location of the apparent horizon (thick dashed line) and 
of the effective excised region (area filled with squares) are 
shown. By this time the star has flattened considerably, all of 
the matter near the rotation axis has fallen inside the appar- 
ent horizon, but a disc of low-density matter remains near the 
equatorial plane, orbiting at very high velocities > 0.2 c. The 
appearance of an effective barrier preventing a purely radial 
infall of matter far from the rotational axis may be the conse- 
quence of the larger initial angular momentum of the this col- 
lapsing matter and of the pressure wave originating from the 
faster collapse along the ritational axis. Note, in fact, that the 
radial velocity at the equator does not increase significantly at 
the stellar surface between t ~ 0.49 and t ~ 0.67 ms, but that 
it actually slightly decreases (cf. the {x, z) planes in the mid 
and lower panels of Fig.|5}- This is the opposite of what hap- 
pens for the radial velocity of the fluid elements in the stellar 
interior on the equatorial plane: it grows also in this time in- 
terval. A more detailed discussion of this behaviour will be 
made in Section fV Dl 

Note that the disc formed outside the apparent horizon is 
not dynamically stable and, in fact, it rapidly accretes onto the 
newly formed black hole. This is shown in Fig.|6l which of- 
fers a magnified view at a later time t = 0.79 ms. At this stage 
the disc is considerably flattened but also has large radial in- 
ward velocities which induce it to be accreted rapidly onto the 
black hole. Note that as the area of the apparent horizon in- 
creases, so does the excised region, which is allowed to grow 
accordingly. This can be appreciated by comparing the areas 
filled with squares in the lower panels of Fig. |5] (referring to 
t = 0.67 ms) with the corresponding ones in Fig. |6l (referring 
to t = 0.79 ms). 

By a time t — 0.85 ms, essentially all (i.e. more than 
99.9%) of the residual stellar matter has fallen within the 
trapped surface of the apparent horizon and the black hole thus 
formed approaches the Kerr solution (see Section |V}. Note 
that a simple kinematic explanation can be given for the insta- 
bility of the disc formed during this oblate collapse and comes 
from relating the position of the outer edge of the disc when 
this first forms, with the location of the ISCO of the newly 
formed Kerr black hole. Measuring accurately the mass and 
spin of the black hole reveals, in fact, that the ISCO is lo- 
cated at a; = 11.08 km, which is always larger than the outer 
edge of the disc (cf. lower panels of Fig.|5jl- Such behaviour 
is not surprising since we are here dealing with initial models 
with a moderate J /NP, that collapse essentially on a dynami- 
cal timescale, and for which pressure gradients cannot play an 
important role. As a result, simple point-like particle motion 
in stationary spacetimes is a sufficient approximation to the 
dynamics. 

Also for model D4, a more detailed discussion (e.g. of the 
evolution of the distribution of angular velocity, or of the disc 
rest mass) will be presented in Section lTV CI Here, however, 
we can anticipate that when analysed more closely the rest 
mass in the disc and outside the apparent horizon is effectively 
very small and that the angular velocity shows appreciable de- 
partures from a uniform profile. 

A more quantitative description of the rest-mass density 
evolution is presented in Fig. |8l where different lines show 
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FIG. 8: Rest-mass density of model D4 normalized to the initial 
value at the stellar centre. The profiles are measured along the x-axis 
on the equatorial plane and refer to different times (see main text for 
details). Line 5, shown as dotted, corresponds to the time when the 
apparent horizon is first found. The inset shows a magnified view of 
the final stages of the evolution using a logarithmic scale and also the 
location of the excised region as it grows in time. 



the profiles of the rest-mass density along the x-axis on the 
equatorial plane. The values are normalized to the initial 
value at the stellar centre, with different labels referring to 
different times and in particular to t = 0.0 (dashed line), 
0.25,0.40,0.49,0.54,0.65,0.67,0.74,0.79 and 0.89 ms, re- 
spectively. Line 5, furthermore, is shown as dotted and refers 
to the time when the apparent horizon is first formed. Af- 
ter this time, the excised region is cut from the computational 
domain as shown in the inset of Fig. |8] which illustrates the 
final stages of the evolution. Note that as the matter falls into 
the black hole, the apparent horizon increases its radius and 
thus the location of the excised region moves outside. This is 
clearly shown in the inset. Note also that the rest-mass density 
does not drop to zero outside the stellar matter but is levelled 
off to the uniform value of the atmosphere, whose rest-mass 
density is seven orders of magnitude smaller than the initial 
central density. 

It should be remarked that such a tenuous atmosphere has 
no dynamical impact and does not produce any increase of the 
mass of the black hole that can be appreciated in our simula- 
tions. With such rest-mass densities, in fact, it would take a 
time ^ IQ^M to produce a net increase of ~ 1% in the black- 
hole mass. Clearly, these systematic errors are well below the 
truncation errors, even at the highest resolutions. 

The simulation ends 2A.t — 0.91 ms ~ 99 M, when the 
rest-mass density is everywhere at the atmosphere level and 
the violations of the Hamiltonian constraint are large. By this 
time the evolution has been carried for more than 28% of the 
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FIG. 9: Evolution of the mass fraction versus time during the collapse of Dl (left panel) and D4 (right panel). The rest mass is measured 
within two-spheres of coordinate radii ri — 0.2, 0.4 and 0.6 times the initial stellar equatorial circumferential radius Re (cf. TableQ}. Marked 
with filled dots on the different lines are the times at which the apparent horizon is first found (the data refers to a simulation with 96^ x 48 
gridzones). The insets in both panels show, on a logarithmic scale, the evolution of the normalized rest mass outside the apparent horizon. 
Note that this is appreciably non-zero for a rather long time in case of model D4. 



total time using a singularity excising region. Also in this 
case, we do not find evidence of shock formation nor of sig- 
nificant deviations from axisymmetry. 



As mentioned in the Introduction, all simulations to-date 
agree that no massive and stable discs form for initial mod- 
els of neutron stars that are uniformly rotating and when a 
polytropic EOS with 1 < < 2 is used. Our results 
corroborate this view and in turn imply that the collapse of 
a rapidly rotating old and cold neutron star cannot lead to 
the formation of the central engine believed to operate in a 
gamma-ray burst, namely a rotating black hole surrounded by 
a centrifugally-supported, self-gravitating torus. Relativistic 
simulations with more appropriate initial data, accounting in 
particular for the extended envelope of the massive progeni- 
tor star which is essential in the so-called collapsar model of 
gamma-ray bursts IstIi . will be necessary to shed light on the 
mechanism responsible for such events. 

Convincing evidence has recently emerged f2§\ that a mas- 
sive disc can be produced if the stellar models are initially 
rotating differentially and with initial total angular momenta 
J/M^ > 1, as it may happen for young and hot neutron stars. 
In this case, the massive disc could emit intense gravitational 
radiation either through its oscillations |88| or as a result of 
the fragmentation produced by non-axisymmetric instabili- 
ties 1 28]. We are presently investigating this possibility and 
the results of our investigation will be reported in a forthcom- 
ing paper 



C. Disc formation and differential rotation 

We now discuss in more detail two interesting properties 
of the matter dynamics in both slowly and rapidly rotating 
models: the evolution of the rest mass outside the apparent 
horizon and the development of differential rotation during 
the collapse. 

In order to monitor the changes of the rest-mass distribu- 
tion during the collapse we define the rest mass within a two- 
sphere of coordinate radius ri < Rf, as (see, for instance, 1891 ) 

M,{n)^ ( pau"V7dV, (4.1) 

J r' KVi 

where d'^x is the 3D coordinate volume element. Shown in 
Fig.|9]is the evolution of the rest masses measured within sev- 
eral representative two-spheres for models Dl (left panel) and 
D4 (right panel), respectively. Different lines refer to differ- 
ent coordinate radii for the two-spheres {i.e. = 0.2, 0.4 
and 0.6 Re, where Re is the initial equatorial circumferen- 
tial radius) and are normalized to the total rest mass within 
the computational domain Af*, shown as a solid line. Marked 
instead with filled dots are the values of M*(ri) at the times 
when the apparent horizon is first found; for simplicity, the 
data shown in Fig.|9]refer to a simulation with 96^ x 48 grid- 
zones, but for this quantities higher resolutions have just the 
effect of shifting the time at which the apparent horizon is first 
found. 

As mentioned before, the excised region is not introduced 
immediately after the apparent horizon has been found, but 
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FIG. 10: Evolution of the averaged angular velocity (Q) (left panel) and of the averaged angular momentum per unit mass (hu^) (right panel). 
Both quantities are measured at the stellar equator, are normalized to the initial value at the stellar surface and refer to both models Dl (upper 
parts) and D4 (lower parts). 



only when this has grown to a sufficiently large size. When 
this happens, the inner part of the computational domain is re- 
moved and the integrals ( 14. 1> are no longer meaningful. As 
a result, all of the curves in Fig.|9lare truncated at the time 
when the excision region is first introduced, which occurs at 
t = 0.72 ms and t — 0.79 ms for models Dl and D4, respec- 
tively. 

A rapid comparison between the two panels of Fig.|9]is suf- 
ficient to identify the differences in the rest-mass evolution 
in slowly and rapidly rotating models. Firstly, the rest-mass 
distribution is very different already initially, being more uni- 
form in Dl and more centrally concentrated in D4, as can 
be appreciated by comparing at — 0.4i?e and 0.6i?e. 
Secondly, the rest-mass infall is much faster for the slowly 
rotating model Dl, while it is more progressive for model 
D4, as shown by the change in the fractional mass ratio at 
Vi — OARe- Finally, the amount of matter outside r.i = 0.4i?e 
at the time when the apparent horizon is found, and which is 
very close to the amount of matter outside the apparent hori- 
zon, is different in the two cases, being essentially zero for 
model Dl and a few percent for model D4. A clearer view 
of this is presented in the two insets of Fig. |5] which show, 
on a logarithmic scale, the evolution of the normalized rest 
masses outside the apparent horizons, i.e. M,(r > r^jj)/M,, 
since these first form and as they grow in time. It is interest- 
ing to note the different behaviour in this case with a rapid 
decrease when the rotation rate is small and a much slower 
one in the case of a rapidly rotating progenitor (Note that the 
two insets cover the same timescale although they refer to a 
different time interval.). 

Two additional comments are worth making. The first one 
is that A/, effectively includes also the rest mass in the atmo- 
sphere but this is always < 10^"* of the total rest mass. The 
second one is that A/* in Fig. |9] does not simply refer to the 



initial value of the total rest mass but is effectively computed 
at each step and appears constant in time because of the ability 
of the code to conserve rest mass. A closer look at the solid 
curve in Fig. |5]re veals, in fact, that A/» varies over time to less 
than one part in 10''. 

An interesting question to ask at this stage is whether these 
uniformly rotating models will develop any degree of differ- 
ential rotation as the collapse proceeds. Part of the interest 
in this comes from the fact that neutron stars are thought to 
rotate differentially, at least during the initial stages of their 
life. This is expected to hold both when the neutron star is 
produced through a stellar core collapse, in which case the 
differential rotation may be present already in the stellar pro- 
genitor and is then amplified during collapse |90], and when 
the neutron star is the end-result of a binary merger of neutron 
stars |91 ]. However, as the neutron star cools and grows older, 
dissipative viscous effects or the coupling with non-turbulent 
magnetic fields are expected to bring the star into uniform ro- 
tation (see 93, 94] for a detailed description of this pro- 
cess). It is therefore interesting to investigate whether a de- 
gree of differential rotation will be produced also during the 
final collapse of a uniformly rotating star to a Kerr black hole. 
To answer this question we have monitored both the averaged 
angular velocity (17), defined as 




(4.2) 



and the corresponding averaged angular momentum per unit 
mass (hucf,), which is a conserved quantity along the path lines 
of fluid elements in an axisymmetric (but not necessarily sta- 
tionary) spacetime 0]. Note that u'^/u* = au"^ — Z?*^ and the 
average over the two different directions is here used to com- 
pensate the small errors that are produced in the evaluation of 
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these quantities near the axes. 

We note that our measure of the differential rotation will 
depend on the specific slicing chosen. However, for the simu- 
lations reported here, the lengthscale of variation of the lapse 
function at any given time is always larger than the stellar ra- 
dius at that time, ensuring that the events on the same timeslice 
are also close in proper time. A useful measure of the differ- 
ential rotation that develops during collapse is the departure 
from unity of the ratio of the values of fl at the centre and at 
the surface of the star on the equatorial plane and it is instruc- 
tive to compare how this varies in the dynamics of the two 
models Dl and D4, which have been evolved using the same 
slicing. 

The time evolution of (f2) and (/iw^) is presented in the two 
panels Fig.^| whose lower parts refer to model D4 while the 
upper ones refer to model Dl. Both quantities are shown nor- 
malized to their initial value at the stellar surface. Let us con- 
centrate on the slowly rotating model first. The different lines 
refer to three representative times which aiet = 0.0 (shown as 
dashed), t = 0.45 and 0.52 ms, respectively. Initially, the an- 
gular velocity is, by construction, uniform throughout the star 
(left panel) and the corresponding specific angular momentum 
grows linearly with the distance from the stellar centre (right 
panel). As the collapse proceeds and the stellar size decreases, 
the angular velocity is expected to increase while the angular 
momentum per unit mass remains constant. This is indeed 
what happens for model Dl, whose specific angular momen- 
tum is conserved with an overall error at the stellar surface 
which is always less than 10% and which decreases with res- 
olution. A similar behaviour is observed also much later in 
the simulation, when the apparent horizon has been found and 
the singularity has been excised. Overall, the angular velocity 
in the collapsing model Dl grows Uke il{t) oc r^T", where 
n ~ 1.5 and therefore less than it would do in the case of 
the collapse of a Newtonian, uniform density star (i.e. n — 2); 
which is a result of relativistic and rotational effects (see ll96ll '). 

A comparison of the lower parts of the two panels in Fig. 1101 
is sufficient to realize that the evolution of the angular veloc- 
ity is rather different for a rapidly rotating stellar model. The 
different lines in this case refer iot — 0.0, 0.48 and 0.65 ms, 
respectively, and it is apparent that a non-negligible degree of 
differential rotation develops as the collapse proceeds, with a 
difference of a factor ^ 2 between the angular velocity of the 
inner and outer parts of the collapsing matter as the apparent 
horizon first appears. Clearly, this differential rotation is pro- 
duced very rapidly and will persist only for a very short time 
before the star is enclosed in a trapped surface. 

It is difficult to establish, at this stage, whether the differ- 
ential rotation generated in this way could produce a phe- 
nomenology observable in some astrophysical context and 
more detailed investigations, in particular of the coupling of 
this differential rotation with magnetic fields fff% losll . are 
necessary. Finally, it is worth remarking that while differen- 
tial rotation develops for model D4 but not for D 1 , the specific 
angular momentum is conserved to the same accuracy in both 
models. 



V. DYNAMICS OF THE HORIZONS 

In order to investigate the formation of a black hole in our 
simulations, we have used horizon finders, available through 
the Cactus framework, which compute both the apparent 
horizon and the event horizon. The apparent horizon, which 
is defined as the outermost closed surface on which all out- 
going photons normal to the surface have zero expansion, is 
calculated at every time step and its location is used to set 
up the excised region inside the horizon. Specific technical 
details about the handling of the excised region for the fields 
are presented in lfl9ll22ll . while a brief discussion of how the 
hydrodynamical excision is performed in Whisky was pre- 
sented in Section HTCl 

In contrast, the event horizon, which is an expanding null 
surface composed of photons which will eventually find them- 
selves trapped, is computed a posteriori, once the simulation 
is finished, by reconstructing the full spacetime from the 3D 
data each simulation produces. In stationary black-hole sys- 
tems, where no mass-energy falls into the black hole, the ap- 
parent and event horizons coincide, but generally (in dynami- 
cal spacetimes) the apparent horizon lies inside the event hori- 
zon. We have here used the fast solver of Thornburg f76ll to 
locate the apparent horizon at every time step, and the level- 
set finder of Diener |99] to locate the event horizon after the 
simulation has been completed and the data produced is post- 
processed. 

In all cases considered, we have found that the event 
horizon rapidly grows to its asymptotic value after forma- 
tion. With a temporal gap of ^ lOM after the formation 
of the event horizon, the apparent horizon is found and then 
it rapidly approaches the event horizon, always remaining 
within it. With the exception of the initial gap, the horizon 
proper areas as extracted from the apparent and event horizon 
are very close (see e.g., Fig.ll6>. 



A. Measuring the Event-Horizon Mass 

We measure the mass of the newly formed black hole to 
estimate the amount of energy that is emitted as gravitational 
radiation during the collapse. In particular, we do a simple 
energy accounting, comparing the mass of the black hole with 
the ADM mass of the spacetime computed by the initial data 
solver on a compactified grid extending to spatial infinity 1 8 1 ] . 
This value is slightly different ( 1 % in the worst case) from the 
one which is instead computed on the finite domain of our 
computational grid at the initial time and after the constraints 
are solved. The difference between the two values can be used 
to define an "error-bar" for our measure of the black-hole mass 
and hence of the energy in gravitational waves {cf. Fig. I15> . 
Two notes are worth making about this error before we go 
on to discuss how the black-hole mass is actually measured. 
Firstly, the difference between the two masses represents the 
truncation error produced by the finite size of the computa- 
tional domain and is conceptually distinct from the truncation 
error introduced by the finite differencing. While the first is 
assumed to be constant in time, the second in general grows 
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FIG. 1 1 : Convergence of the measure of the black-hole mass as the 
resolution is increased. The curves refer to estimates using the event- 
horizon equatorial circumference [i.e. eq. <5.1H and have been ob- 
tained using 288^ X 144 and 192^ x 96 zones, respectively. Shown 
in the small inset are the results for model Dl, while those for model 
D4 are in the main panel. 
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FIG. 12: Evolution of the event-horizon mass M — Ceq/iir for 
models Dl and D4. Different lines refer to the different initial 
guesses for the null surface and are numbered 0, 1 and 2. Note that 
they converge exponentially to the correct event-horizon surface for 
decreasing times. 



with time (especially after the excision is made) and is moni- 
tored through the calculation of the constraint equations. Sec- 
ondly, this error-bar sets a global lower limit on the accuracy 
of our measure of asymptotic quantities and therefore on the 
energy lost to gravitational waves during the collapse. No reli- 
able measure of this lost energy can be made below the error- 
bar even if the constraint equations are solved to a larger pre- 
cision (this will be discussed in more detail in Section fVCt . 

The first and simplest method of approximating the black- 
hole mass is to note that, for a Kerr (or Schwarzschild) black 
hole, the mass can be found directly in terms of the event- 
horizon geometry as 

M = ^, (5.1) 

where Ceq = J^^ y^g^dcf) is the proper equatorial circumfer- 
ence. Provided there is a natural choice of equatorial plane, 
it is expected that, as the black hole settles down to Kerr, Ceq 
will tend to the correct value. However, as numerical errors 
build up at late times it may be impossible to reach this asymp- 
totic regime with sufficient accuracy. 

The estimate of M coming from the use of i5.H is pre- 
sented in Fig. ^2 which shows the time evolution of the event- 
horizon equatorial circumference. The two lines refer to two 
different resolutions (288'^ x 144 and 192^ x 96 zones, respec- 
tively) and should be compared with the value of the ADM 
mass A/adm (indicated with a short-long dashed line), and 
with the error-bars as deduced from the initial data. Shown in 



the small inset are the results for model Dl, while those for 
model D4 are in the main panel. 

Note that if a measure of the event horizon is not available, 
eq. i5.l\ could be computed using the equatorial circumfer- 
ence of the apparent horizon (this is what was done, for in- 
stance, in 1231 '). Doing so would yield results that are similar 
to those shown in Fig. ^3 although with a slightly larger de- 
viation from A/adm- This is because we have found the ap- 
parent horizon to systematically underestimate the equatorial 
circumference. In particular, in the high-resolution run for 
model D4, the differences between the apparent and event- 
horizon equatorial circumferences are < 2 %. 

Clearly, as the equatorial circumference grows, the agree- 
ment with the expected ADM mass improves as it does with 
the use of higher spatial resolution. However, equally evident 
is that the eiTors grow as the collapse proceeds and this is due, 
in part, to the loss of strict second-order convergence at later 
times, but also to the way the event horizon is found. The 
level-set approach of [99], in fact, needs initial guesses for 
the null surface, which converge exponentially to the correct 
event-horizon surface for decreasing times, hence introduces 
a systematic error in the calculation of the event horizon at late 
times. This is shown in Fig.^] which presents the evolution 
of the event-horizon mass M = Ceg/47r for models Dl and 
D4. Different lines refer to the different initial guesses and 
are numbered "0", "1" and "2", respectively (note that for the 
curves shown in Fig.^Jthe initial guesses "0" and "1" were 
used for cases D4 and Dl, respectively). 
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FIG. 13: Fitting the oblateness of the event horizon to QNMs of a Kerr black hole. The fit is shown with the solid line, while the open circles 
represent the computed values of Cr- The estimate for Cr of a Kerr black hole having the fitted value of a/Mhor is shown with a dashed line. 



B. Measuring the angular momentum of the black hole 

A major difficulty in an accurate measurement of M lies in 
the calculation of its non-irreducible part, i.e. in the part that 
is proportional to the black-hole angular momentum J. We 
now discuss a number of different ways to calculate J from 
the present simulations; these measurements will then be used 
to obtain alternative estimates of M in Section lVCl 



1. Measuring J from tlie Horizon Distortion 

In a series of papers studying the dynamics of distorted 
black-hole spacetimes, it was shown that the horizon geom- 
etry provides a useful measure of th e black-hole properties 
both in vacuum fl0(]l ll0lL ll02'.'l03'l and when these are ac- 
creting matter axisymmetrically 1 104]. In particular, the idea 
is to look at the distortion of the horizon using the ratio of 
polar and equatorial proper circumferences, Cr = Cpoi/Ceq. 
For a perturbed Kerr black hole this is expected to oscillate 
around the asymptotic Kerr value with the form of a quasi- 
normal mode (QNM). By fitting to this mode we extract an 
estimate of the a ngular momentum parameter a/Mhor from 
the relation llOSll 

= v^l- (-1.55 + 2.55a)' , (5.2) 

Mhor 

where we have indicated with Mf^or the black-hole mass as 
measured from expression i5.2\ . which coincides with M only 
if the spacetime has become axisymmetric and stationary. The 
fit through expression ( 15. 2t is expected to be accurate to ^ 
2.5% 1105]. 

The fit itself depends on an initial guess for a/Mhor and 
we start from a Schwarzschild black hole and iterate till the 



desired convergence is reached. This measure is not gauge 
invariant, although eq. ( 15. 2> is independent of the spatial co- 
ordinates up to the definition of the circumferential planes, 
but works adequately with the gauges used here. The fit is 
best performed shortly after black-hole formation as the oscil- 
lations are rapidly damped. This minimizes numerical errors 
but in those cases where matter continues to be accreted, it 
may lead to inaccurate estimates of the angular momentum. 

Examples of the fitting procedure are shown in Fig. [O] 
in which the fit is shown as a solid line, while the open cir- 
cles represent the computed values of Cr', these are slightly 
noisy as a result of the interpolation needed by the level-set 
approach to find points on the horizon two-surface ]99]. The 
estimate for Cr of a Kerr black hole having the fitted value 
of a/Mhor is shown as a dashed line. Note that the values of 
a/Mhor — 0.21 and a/Mhor — 0.54 are very close to the total 
J/M^ of the initial stellar models, i.e. 0.2064 and 0.5433, as 
shown in Tablellll This demonstrates that, to within numerical 
accuracy, the entire angular momentum of the spacetime ends 
up in the black hole. 

Using expression ( 15. 2> to estimate the angular momentum 
J introduces an error, if the black hole has not yet settled to a 
Kerr solution. Having this in mind, however, it is possible to 
estimate the angular momentum as 

2. Measuring J with the dynamical-horizon framework 

A second method of approximating J and hence measuring 
M is to use the isolated and dynamical -honzon fram eworks 
of Ashtekar and collaborators 1106, 107, 108, 109, 110]. This 
assumes the existence of an axisymmetric Killing vector field 
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FIG. 14: Comparison of the different measures of the angular momentum for the cases Dl (left panel) and D4 (right panel). The estimate 
using the fit to the circumference ratio (see left panel of Fig. ll3> is also shown. The dynamical-horizon spin measurement is considerably more 
accurate at late times as the event-horizon surfaces will diverge exponentially at this point. Shown with the horizontal short-long dashed lines 
are the values of (J/M^ )adm in the two cases as measured from the initial data (see main text for details). 



TABLE IT: Estimates of the black-hole angular momentum J /h'P 
through the oblateness of the event horizon. The oscillations in the 
oblateness of the event horizon can be fitted to the normal modes 
of a Kerr black hole. Note that for each model the measured an- 
gular momentum is remarkably close to that of the initial spacetime 
( J/i\f ^)adm. Also reported are the initial ADM mass, the value of 
the equatorial circumference as obtained through the fit (Cr)EH, and 
the corresponding value obtained through the estimated spin param- 
eter (C,.)Koir. 



Model AfADM (J/A'/)iDM 


(J/M2)eh 


(Cr)EH 


(Cr) Kerr 


Dl 1.6653 0.2064 


0.21 


0.99 


0.9916 


D2 1.7281 0.3625 


0.36 


0.97 


0.9734 


D3 1.7966 0.4685 


0.47 


0.95 


0.9544 


D4 1.8606 0.5433 


0.54 


0.94 


0.9372 



intrinsic to a marginally trapped surface such as an apparent 
horizon. This gives an unambiguous definition of the spin 
of the black hole and hence of its total mass. If there is an 
energy flux across the horizon, the isolated-horizon frame- 
wor k needs to be extended to the dynamica l-hohzon formal- 
ism j[^[lTl. 

In practice, our approach to the dynamical-horizon frame- 
work has been through the use of a code by Schnetter which 
implements the algorithm of fl lO*] to calculate the horizon 
quantities. The advantage of the dynamical-horizon frame- 
work is that it gives a measure of mass and angular momentum 
which is accurately computed locally, without a global recon- 
struction of the spacetime. One possible disadvantage is that 
the horizon itself is required to be (close to) axisymmetric; so 
in case it deviates largely from axial symmetry, no informa- 



tion can be found. However, because arbitrarily large distor- 
tions are allowed as long as they are axisymmetric, we have 
not encountered problems in applying the dynamical-horizon 
framework to the horizons found in our simulations. 



3. Measuring J from the Angular Velocity of the Event Horizon 

A third method for computing J only applies if an event 
horizon is found and if its angular velocity has been measured. 
In a Kerr background, in fact, the generators of the event hori- 
zon rotate with a constant angular velocity ut = ~gtct>/gct>ct> = 
\/gtt/ g4>ct>- In this case the generator velocity can be directly 
related to the angular momentum parameter as 



a 

M 



J 

M2 



47r 



1/2 



(5.4) 



As with the previous approximations, expression \5A\ is 
strictly valid only for a Kerr black hole and will therefore con- 
tain a systematic error which, however, decays rapidly as the 
black-hole perturbations are damped. On the other hand, the 
event horizon generator velocities have the considerable ad- 
vantage that everything is measured instantaneously and the 
values of u are valid whether or not the background is an iso- 
lated Kerr black hole. The disadvantage, though, is that, as 
mentioned above, the numerical event horizon surfaces be- 
come systematically less accurate at late times (cf. Fig.ll2>. 
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Comparison of angular momentum measurements 



and obtain 



A detailed comparison of the three different methods for 
measuring the angular momentum of the black hole is shown 
in Fig. [2] The measurement of angular momentum using 
the angular velocity of the generators is shown as solid lines. 
Both for slowly (left panel) and rapidly (right panel) rotating 
stellar models, the event horizon has zero area (and thus zero 
angular momentum) when it is first formed. However, as the 
rotating matter collapses, the event horizon-area and angular 
momentum grow, the black hole is spun up and, to numeri- 
cal accuracy, the total angular momentum of the spacetime is 
contained within the black hole {cf. Fig.ll3>. At late times, the 
estimate using the generator velocities of the event horizon 
drifts away, probably due to a combination of gauge effects 
and the systematic errors in the trial guesses for the null sur- 
faces. 

In the case of the slowly rotating model Dl, in particular, 
the estimate from the dynamical-horizon finder is perfectly 
stable {cf. dashed line in the left panel of Fig. I14t . indicat- 
ing that an approximately stationary Kerr black hole has been 
formed by the time the simulation is terminated. In the case 
of the rapidly rotating model D4, however, this is no longer 
the case as matter continues to be accreted also at later times, 
when the errors have also increased considerably. As a result, 
the measure of the spin through the dynamical-horizon finder 
is less accurate and does not seem to have stabilized by the 
time the simulation ends (cf. dashed line in the right panel of 
Fig. I14> . This may indicate that the final black hole has not 
settled down to a Kerr black hole on the timescales considered 
here. 



C. black-hole mass from the Christodoulou formula 

It was shown by Christodoulou that, in the axisymmetric 
and stationary spacetime of a Kerr black hole, the square of 
the black-hole mass A/ is given by 



(5.7) 
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where is the irreducible mass, A is the event-horizon 
proper area, and J is the black-hole angular momentum. As 
the black hole approaches a stationary state at late times, the 
apparent and event horizons will tend to coincide and in that 
case the mass of the black hole is very well approximated by 
the above formula. 

We have applied the above formula, using the various meth- 
ods for measuring the angular momentum J. In particular, us- 
ing the method for obtaining J from the distortion of the event 
horizon, through eq. \5.3l . the black-hole mass is given by 
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If, on the other hand, J is found from the angular velocity oj 
of the event horizon, then it is possible to use \5.A\ in \5.5\ 
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In the framework of dynamical horizons, expression \5.5\ 
holds for any axisymmetric isolated or dynamical horizon, in- 
dependently of whether it is stationary. 

Figure^]collects the four different ways of measuring the 
black-hole mass for the collapse of models Dl and D4. The 
different lines refer to the different approaches we have out- 
lined above. In addition, we display the irreducible mass Afirr. 
The left panel of Fig. ^] in particular, shows the results of 
the different measures for the slowly rotating model Dl. Be- 
cause in this case all of the matter rapidly collapses into the 
black hole, the different estimates of the total mass agree very 
well. However, already in this slowly rotating case the irre- 
ducible mass of the apparent horizon is noticeably lower The 
left panel also shows that while the different methods provide 
comparable estimates, only the one corresponding to eq. ( 15. Il l 
(i.e. the solid line) falls for some time within the error-bar 
provided by the initial estimate of A/adm (this is particularly 
evident in the inset). Because when this happens the norms 
of the Hamiltonian constraint have not yet started to grow ex- 
ponentially and the largest value of the constraint violation is 
about an order of magnitude smaller (i.e. the Li norm of the 
Hamiltonian constraint is ~ 4.9 x 10""* alt — 0.56 ms) we 
can use the error-bar in A/adm to place an upper bound of 
0.5% A/adm to the energy lost through the emission of grav- 
itational radiation in this case. Clearly, the true bound is cer- 
tainly considerably lower and we expect that with accuracies 
comparable to the ones of 2D simulations, our estimates of the 
efficiency of gravitational radiation emission could converge 
to the values of Stark and Piran 1 32|. 

The right panel of Fig. on the other hand, shows the 
results of the different mass measures for the rapidly rotat- 
ing model D4. In this case, the contribution from the spin 
energy is considerably larger and noticeable differences ap- 
pear among the different approaches. Since all seem to have 
systematic errors, this makes it less trivial to establish which 
method is to prefer. On one hand, those methods using in- 
formation from the event-horizon equatorial circumference or 
that fit the perturbations of the event horizon [i.e. eqs. ( 15.11 1 
and ( 15. 6H seem to provide accurate estimates at earlier times 
but suffer of the overall inaccuracy at later stages, when the 
initial guesses for the null surface are distinct. It is indeed 
at these early times that these measurements are within the 
error-bar provided by the initial estimate of A/adm. On the 
other hand, those methods that measure the angular velocity 
of the null generators [i.e. eq. \5.1V \ or that use the dynamical- 
horizon framework, produce reasonably accurate estimates, 
that converge with resolution, that monotonically grow in time 
and that are within the error-bar of the initial estimate of 
.^^ADM. Furthermore, in the case of the dynamical-horizon 
framework, this is not only physically expected, given that 
a small but non-zero fraction of the matter continues to be 
accreted nearly until the end of the simulation, but it is also 
guaranteed analytically. 

Because of these differences in the measures of A/ and be- 
cause the black hole does not have time to settle down to a 
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FIG. 15: Comparison of the five different approaches used in the measure of the total black-hole mass for the collapse of models Dl and D4. 
Different lines refer to the different methods discussed in the main text. The left panel (model Dl), shows that the different methods are overall 
comparable when the rotation is slow, but that differences are already present (this is as shown in the inset). The right panel (model D4) shows 
that the different measures can be considerably different when the rotation is large. 



constant total mass, the upper bound on the energy emission is 
more conservative than in case Dl. In particular, taking again 
as a reference the time when the estimate relative to eq. ( 15. 1> 
is within the error-bar (i.e. at t ~ 0.70 ms) and the largest 
value of the constraint violation is about an order of magni- 
tude smaller (i.e. the Li norm of the Hamiltonian constraint 
is ~ 1.2 X 10^"^) and is not yet growing exponentially, we 
place an upper bound of 1% Madm on the energy lost through 
gravitational radiation. Once again, we expect the true value 
to be considerably smaller. 

One obvious and expected result is that the irreducible mass 
in the collapse of model D4 (the dot-dashed line in the right 
panel of Fig.ll5ll deviates by a large amount from the actual 
black hole mass, since it does not include the rotational energy 
of the black hole. 

Finally, we will make a comment on the different methods 
used for measuring the mass and spin of a black hole in a nu- 
merical simulation. Although the direct comparison of many 
different methods employed here have provided valuable in- 
formation on the dynamics of the system, we have found the 
dynamical-horizon framework to be simple to implement, ac- 
curate and not particularly affected by the errors from which 
equivalent approaches seem to suffer, as shown in our Figs. 1151 
and ^5 As a result, we recommend its use as a standard tool 
in numerical relativity simulations. 



D. Reconstructing the global spacetime 

All of the results presented and discussed in the previous 
Sections describe only a small portion of the spacetime which 
has been solved during the collapse. In addition to this, it is 
interesting and instructive to collect all of these pieces of in- 



formation into a global description of the spacetime and look 
for those features which mark the difference between the col- 
lapse of slowly and rapidly rotating stellar models. As we dis- 
cuss below, these features emerge in a very transparent way 
within a global view of the spacetime. 

To construct this view, we use the worldlines of the most 
representative surfaces during the collapse, namely those of 
the equatorial stellar surface, of the apparent horizon and of 
the event horizon. For all of them we need to use properly de- 
fined quantities and, in particular, circumferential radii. The 
results of this spacetime reconstruction are shown in Fig. 1161 
whose left and right panels refer to the collapse of models Dl 
and D4, respectively. The different lines indicate the world- 
lines of the circumferential radius of the stellar surface (dotted 
line), as well as of the apparent horizon (dashed line) and of 
the event horizon (solid line). Note that for the horizons we 
show both the equatorial and the polar circumferential radii, 
with the latter being always smaller than the former. For the 
stellar surface, on the other hand, we show the equatorial cir- 
cumferential radius only. This is because the calculation of 
the stellar polar circumferential radius requires a line integral 
along the stellar surface on a given polar slice. Along this 
contour one must use a line element which is suitably fitted to 
the stellar surface and diagonalized (see 1 103] for a detailed 
discussion). In the case of model D4, however, this is difficult 
to compute at late times, when the disc is formed and the Une 
integral becomes inaccurate. 

Note that in both panels of Fig.^]the event horizon grows 
from an essentially zero size to its asymptotic value. In con- 
trast, the apparent horizon grows from an initially non-zero 
size and, as it should, is always contained within the event 
horizon. At late times, the worldlines merge to the precision 
at which we can compute them. A rapid look at the two pan- 
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FIG. 16: Evolution of the most relevant surfaces during the collapse for the cases Dl and D4. Solid, dashed and dotted lines represent the 
worldlines of the circumferential radii of the event horizon, of the apparent horizon and of the stellar surface, respectively. Note that for the 
horizons we plot both the equatorial and the polar circumferential radii, while only the equatorial circumferential radius is shown for the stellar 
surface. Shown in the insets are the magnified views of the worldlines during the final stages of the collapse. 



els of Fig.^]is sufficient to appreciate the different properties 
in the dynamics of the collapse of slowly and rapidly rotating 
models. 

Firstly, in the case of model Dl, the differences between the 
equatorial and polar circumferential radii of the two trapped 
surfaces are very small and emerge only in the inset that offers 
a magnified view of the worldlines during the final stages of 
the collapse. This is not the case for model D4, for which 
the differences are much more evident and can be appreciated 
also in the main panel. Of course, this is what one expects 
given that the ratio of these two quantities depends on a/M 
and is ~ 1 for a slowly rotating black hole {cf. TableHH. 

Secondly, the worldlines of the stellar equatorial circumfer- 
ential radius are very different in the two cases. In the slowly 
rotating model Dl, in particular, the star collapses smoothly 
and the worldline always has negative slope, thus reaching 
progressively smaller radii as the evolution proceeds (cf. left 
panel of Fig. [T6l. By time t ~ 0.59 ms, the stellar equato- 
rial circumferential radius has shrunk below the correspond- 
ing value of the event horizon. In the case of the rapidly ro- 
tating model D4, on the other hand, this is no longer true and 
after an initial phase which is similar to the one described for 
Dl, the worldline does not reach smaller radii. Rather, the 
stellar surface slows its inward motion and, at around t ^ 0.6 
ms, the stellar equatorial circumferential radius does not vary 
appreciably. Indeed, the right panel of Fig. [16| shows that at 
this stage the stellar surface moves to slightly larger radii. This 
behaviour marks the phase in which a flattened configuration 
has been produced and the material at the outer edge of the 
disc experiences a stall {cf. the middle and lower panels of 
Fig.|5jl. As the collapse proceeds, however, also this material 
will not be able to sustain its orbital motion and, after t ^ 0.7 
ms, the worldline moves to smaller radii again. By a time 



t ~ 0.9 ms, the stellar equatorial circumferential radius has 
shrunk below the corresponding value of the event horizon. 



VI. CONCLUSION 

Although 3D numerical relativity has been a very active re- 
search area for several years now, there are still a number of 
technical issues to be addressed and physical problems to be 
investigated in detail. Separate progress has been made so 
far in obtaining long-term stable evolutions of vacuum space- 
times and of spacetimes with matter. Both of them have posed 
significant numerical problems because of the existence of 
horizons containing physical singularities, in one case, and the 
development of non-linear hydrodynamical phenomena such 
as shocks, in the other In black-hole vacuum spacetimes, 
these problems have successfully been dealt with by using bet- 
ter suited formulations of the Einstein equations and by em- 
ploying excision techniques for the regions of the spacetime 
containing the singularity. In spacetimes containing matter, 
on the other hand, sophisticated numerical techniques (such as 
the HRSC methods) have been employed to accurately track 
the dynamics of the shocks. 

Here, we have combined these two different approaches 
by implementing excision techniques within a forming hori- 
zon, thus following the dynamics of the matter as it accretes 
onto the developing black hole. We have shown that doing so 
allows the numerical evolution to proceed uninhibited from 
fully regular initial conditions of matter in equilibrium and 
devoid of trapped surfaces, up to a vacuum spacetime featur- 
ing an event horizon enclosing an excised physical singularity. 
This new important ability in numerical relativity evolutions 
will help in a more detailed investigation of complex astro- 
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physical systems, such as the coalescence of neutron star bi- 
naries, considered as a prime candidate for the detection of 
gravitational waves, and of the collapse of stellar cores, con- 
sidered as the progenitors of gamma-ray bursts. 

Our hydrodynamical excision technique is implemented 
within a new 3D general-relativistic numerical evolution 
code that combines state-of-the-art numerical methods for the 
spacetime evolution {i.e. the NOK formulation of the Einstein 
equations with Gamma-driver shift conditions) with an accu- 
rate hydrodynamical evolution employing several high-order 
HRSC methods. The evolution of the spacetime and of the 
hydrodynamics is coupled transparently through the method 
of lines, which allows for the straightforward implementation 
of various different time-integrators. 

As a first astrophysical problem for this novel setup, we 
have here focused on the collapse of rapidly rotating relativis- 
tic stars to Kerr black holes. The stars are assumed to be in 
uniform rotation and dynamically unstable to axisymmetric 
perturbations. While the collapse of slowly rotating initial 
models proceeds with the matter remaining nearly uniformly 
rotating, the dynamics is shown to be very different in the case 
of initial models rotating near the mass-shedding limit, for 
which strong differential rotation develops. Although the stars 
become highly flattened during collapse, attaining a disc-like 
shape, the collapse cannot be halted because the specific an- 
gular momentum is not sufficient for a stable disc to form. In- 
stead, the matter in the disc spirals towards the black hole and 
angular momentum is transferred inward to produce a spin- 
ning black hole. 

Several different approaches have been employed to com- 
pute the mass and angular momentum of the newly formed 
Kerr black hole. Besides more traditional methods involving 
the measure of the geometrical properties of the apparent and 
event horizons, we have fitted the oscillations of the perturbed 
Kerr black hole to specific quasi-normal modes obtained by 
linear perturbation theory. In addition, we have also consid- 
ered the recently proposed isolated and dynamical horizon 
frameworks, finding it to be simple to implement and yield- 
ing estimates which are accurate and more robust than those 
of other methods. This variety of approaches has allowed for 
the determination of both the mass and angular momentum of 
the black hole with an accuracy unprecedented for a 3D sim- 
ulation. These measures, in turn, have allowed us to set upper 
limits on the energy and angular momentum that could be lost 
during the collapse in the form of gravitational radiation. 

Work using mesh-refinement techniques is already in 
progress to extract more precise information and waveforms 
for the gravitational radiation recorded at large distances from 
the collapsing stars. This aspect of the gravitational collapse 
has not yet been considered in full 3D simulations and will be 
reported in a forthcoming paper f^S*]. Finally, all of the tech- 
niques discussed here will also be applied to the study of the 
collapse of differentially rotating stars, governed by more re- 
alistic and non-isentropic EOSs. The expectation is that initial 
data with J/AP > 1 can be constructed in this case, whose 
collapse could lead to the formation of a massive disc orbiting 
around the newly formed Kerr black hole L12..28»..38il . 



APPENDIX A: NUMERICAL METHODS 

In this Appendix we focus on the numerical methods that 
the Whisky code incorporates for the solution of the gen- 
eral relativistic hydrodynamics equations. The corresponding 
methods for the spacetime equations are those implemented in 
the Cactus code and they have been reported elsewhere and 
the interested reader is addressed to LI 3.. .57.1 for more details. 

As mentioned in the main text, our code uses high- 
resolution shock-capturing methods based on reconstruction 
evolution methods. In this approach, all variables q are repre- 
sented on the numerical grid by cell-integral averages. The 
function is then reconstructed within each cell, usually by 
piecewise polynomials in a way that preserves conservation 
of the variables q. This gives two values at each cell bound- 
ary which are then used to solve (approximately) the Riemann 
problem, giving the flux through the cell boundary. A Method 
of Lines approach is then used to update in time. We will here 
give brief descriptions of each method, but further details can 
be found in [39]. 



1. Time update: the method of lines 

The reconstruction methods guarantee that a prescribed or- 
der of accuracy is retained in space. However, the need to 
retain a high-order accuracy also in time can complicate con- 
siderably the evolution from a time-level to the following one. 
As a way to handle this efficiently, we have chosen to fol- 
low a MoL approach fTS','??]. Here, the continuum equations 
are considered to be discretized in space only. The resulting 
system of ordinary differential equations (ODEs) can then be 
solved numerically with any stable solver. This method mini- 
mizes the coupling between the spacetime and hydrodynamics 
solvers and allows for a transparent implementation of differ- 
ent evolution schemes. 

MoL itself does not have a precise truncation error but, 
rather, it acquires the truncation order of the time-integrator 
employed. Several integrators are available in our imple- 
mentation of MoL, including the second-order Iterative Crank 
Nicholson (ICN) solver and Runge-Kutta (RK) solvers of first 
to fourth-order accuracy. The second and third-order RK 
solvers are known t o be T VD whilst the fourth-order is known 
to not be TWD\V&Jl T?]. As the coupling between the space- 
time and the hydrodynamics is only second-order accurate, we 
typically use the ICN solver 

The calculation of the right hand side to feed to the ODE 
splits into the following parts: 

L Calculation of the source terms s(q(x^^'' , a;*;^' a; l'^'')) at 
all the grid points. 

2. For each direction a:*^*-' : 

• Reconstruction of the data q to both sides of a cell 
boundary. In this way, two values q^^ and q^ of 
qj.+i/2 are determined at the cell boundary; q^ 
is obtained from cell ji (left cell) and q^ from 
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cell ji + 1 (right cell) (see Appendix IA 2l for more 
details). 

• Solution at cell boundary of the approximate Rie- 
mann problem having the values ^ as initial 
data (see Appendix IA 3l for more details). 

• Calculation of the inter-cell flux f^^' (iji+i/a); 
that is, of the flux across the boundary between a 
cell (e.g., the ji-th) and its closest neighbour (e.g., 
the (j, + l)-th). 

3. Recovery of the primitive variables and computation of 
the stress-energy tensor for use in the Einstein equa- 
tions. 

As a result of steps 1.-3., the core of the Whisky code 
is effectively represented by two routines. One that recon- 
structs the function q at the boundaries of a computational 
cell and another one that calculates the inter-cell flux f at this 
cell boundary. 



2. Reconstruction methods 

For the reconstruction procedure we have implemented sev- 
eral different approaches, including slope-limited TVD meth- 
ods, the Piecewise Parabolic Method |48] and the Essentially 
Non-Oscillatory method \49]. 

The TVD method uses limiters to avoid oscillations at 
shocks: we typically use the Van Leer monotized centred 
method, although a variety of others (minmod, Superbee) 
have also been implemented. This method is simple and 
computationally the least expensive to implement, but is at 
most second-order accurate, dropping to first-order at local 
extrema. 

The PPM method of Colella and Woodward il is a com- 
posite reconstruction method that has special treatments for 
shocks, where the reconstruction is modified to retain mono- 
tonicity, and contact surfaces, where the reconstruction is 
modified to sharpen the jump. PPM contains a number of tun- 
able parameters, but we always use those suggested by 14^ . 
PPM is third-order accurate at most. 

The ENO methods have a large number of variants, with 
the common theme that the "least oscillatory" stencil amongst 
all possible stencils of a given order is chosen. In practice we 
use it in its simplest form: direct piecewise polynomial recon- 
struction of the variables, as described in jl 15ll . ENO has no 
tunable parameters besides the order of accuracy, which may 
be arbitrary. 

All these methods are stable in the presence of shocks. By 
default we use PPM as this seems to be the best balance be- 
tween accuracy and computational efficiency, as shown, for 
example, in II84I1 . and is our best choice for all the test evolu- 
tions we present in AppendixlHI 



3. Approximate Riemann problem solvers 

Once a reconstruction procedure has provided data on ei- 
ther side of each cell boundary, this is then used to specify 
the initial states of the semi-infinite piecewise constant Rie- 
mann problems. Since the exact solution of the Riemann 
problem 1 1 16j| is still too costly to use, even when recast in 
an efficient form lllTllllSil . we have here implemented three 
different approximate Riemann solvers. 

The first and simplest method implemented is the HLLE 
method [501 . This approximates the solution by only two 
waves with the intermediate state given by the conservation 
of the mass-flux. This method is very efficient but diffusive. 
The second method is the Roe solver 1 119|. This solves a lin- 
earized problem at each boundary, approximating every wave 
by either a shock or a contact discontinuity. This method is 
less efficient but very accurate. However, it may have prob- 
lems near sonic points. The third method is the Marquina 
solver H20i il2L] . This is similar to the Roe solver, except 
that at possible sonic points a Lax-Friedrichs flux (analogous 
to the HLLE method) is used, ensuring that the solution does 
not contain rarefaction shocks. Note that we use the modified 
method of lIsTIl instead of the original method. 

Both the Roe and Marquina solvers require the computation 
of the eigenvalues and eigenvectors (from both the right and 
left cell) of the linearized Jacobian matrices A^^ and given 
by = -A-j^q^ and = A^q^. We use an implemention 
of the analytic expression for the left eigenvectors |52], thus 
avoiding the computationally expensive inversion of the three 
5x5 matrices of the right eigenvectors, associated to each spa- 
tial direction. We also use a compact version of the flux for- 
mula (a variant on the methods described in |53]) to increase 
speed and accuracy. These improvements bring a ~ 40% 
reduction of the computational time spent in the solution of 
the hydrodynamics equations and a ~ 5 — 10% reduction in 
evolutions involving also the time integration of the Einstein 
equations. The small overall gain in efficiency is due to the 
fact that only around half of the computational time is spent 
computing the update of the hydrodynamic variables, with the 
other half spent in the update of the spacetime field variables. 
Of the time spent computing the hydrodynamic update terms 
around one third is spent solving the Riemann problem. These 
improvements were made by Joachim Frieben at the Univer- 
sidad de Valencia. 



4. Treatment of the atmosphere 

At least mathematically, the region outside our initial stel- 
lar models is assumed to be perfect vacuum. Independently 
of whether this represents a physically reaUstic description 
of a compact star, the vacuum represents a singular limit of 
the equations (j^TOJ and must be treated artificially. We have 
here followed a standard approach in computational fluid- 
dynamics and added a tenuous "atmosphere" filling the com- 
putational domain outside the star. This approach, which was 
implemented also in 1 14,^1'], was instead not used in the sim- 
ulations presented in I.122il . where new strategies were sug- 
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gested. Unfortunately, these corrections would not have been 
effective here because of the conservation form of our hydro- 
dynamical equations ( I2.1Q> . which we have used to guarantee 
the correct evolution of shocks. 

We have treated the atmosphere as a perfect fluid governed 
by the same polytropic EOS used for the bulk matter, but hav- 
ing a zero coordinate velocity. Furthermore, the rest-mass 
density is several (usually 7) orders of magnitude smaller than 
the initial central density. Note that the atmosphere used for 
the calculation of the initial data and the one evolved during 
the simulations need not be the same. Indeed, for the ini- 
tial stellar models used for the collapse calculations we have 
typically set the atmosphere to be two orders of magnitude 
smaller than the evolved one to minimise spurious matter ac- 
cretion onto the black hole. In the pulsation tests presented in 
AppendixlHl on the other hand, the initial and evolved atmo- 
spheres are the same. 

The evolution of the hydrodynamical equations in grid- 
zones where the atmosphere is present is the same as the one 
used in the bulk of the flow. Furthermore, when the rest mass 
in a gridzone falls below the threshold set by the atmosphere, 
that gridzone is simply not updated in time. 

As mentioned in Section lTVBI the use of a tenuous atmo- 
sphere has no dynamical impact and does not produce any 
increase of the mass of the black hole that can be appreciated 
in our simulations. With the rest-mass densities used here for 
the atmosphere, in fact, and using the mass accretion rates 
measured once the apparent surface is first found, we have es- 
timated that a net increase of ^ 1% in the black-hole mass 
would require an integration time of ^ 10"' Af. These system- 
atic errors are well below our truncation errors, even at the 
highest resolutions we can afford. 



APPENDIX B: NUMERICAL TESTS 

Several tests have been performed to assess the stability and 
accuracy of our code, some of which have also provided new 
interesting physical results. Here, we briefly discuss some of 
the most representative of these tests and results, referring the 
interested reader to L39i] for more details. 

First of all, we consider a standard shock-tube test, setting 
as initial data a global Riemann problem, in particular one 
in which the initial discontinuity is orthogonal to the main 
diagonal of our cubic grid. More precisely the initial data 
consists of a left (L) and right (R) states with values given 
respectively by 

PB =1; = 1-666 X 10"^; VR^O , 
PL =10; PL = 13.333; vl = ■ 

In Fig. [n] we show the solution at a given time together with 
the exact solution. The excellent agreement of the two sets of 
curves is particularly remarkable if one bears in mind that the 
evolution is fully 3D and not simply along a coordinate axis. 

Next, we consider the evolution of a stable relativistic poly- 
tropic spherical (TOV) star As this is a static solution, no 
evolution is expected. Yet as shown in Fig. both a small 
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FIG. 17: Solution of a Riemann problem set on the main diagonal of 
the cubic grid. The figure shows the comparison of the hydrodynam- 
ical variables evolved by Whisky, indicated with symbols, with the 
exact solution. The numerical simulation was obtained with the van 
Leer reconstruction method and the Roe solver, on a grid with 140^ 
points. 

periodic oscillation and a small secular increase of the cen- 
tral density of the star are detected during the numerical evo- 
lution of the equations. Both effects have, however, a single 
explanation. Since the initial data contains also a small trunca- 
tion error, this is responsible for triggering radial oscillations 
which appear as periodic variations in the central density. As 
the resolution is increased, the truncation error is reduced and 
so is the amplitude of the oscillation. The secular growth, 
on the other hand, is a purely numerical problem, probably 
related to the violation of the constraint equations. As for 
the oscillation amplitude, also the secular growth converges 
to zero with increasing resolution. Note that the convergence 
rate is not exactly second-order but slightly smaller, because 
the reconstruction schemes are only first-order accurate at lo- 
cal extrema {i.e. the centre and the surface of the star) thus 
increasing the overall truncation eiTor llssll . 

In order to further investigate the accuracy of our imple- 
mentation of the hydrodynamics equations, we have sup- 
pressed the spacetime evolution and solved just the hydrody- 
namics equations in the fixed spacetime of the initial TOV 
solution. This approximation is refeiTed to as the "Cowling 
approximation" and is widely used in perturbative studies of 
oscillating stars. In this case, in addition to the confirmation 
of the convergence rate already checked in fully evolved runs, 
we have also compared the frequency spectrum of the numer- 
ically induced oscillation with the results obtained by an inde- 
pendent 2D code |84| and with perturbative analyses. 

In Fig. ^] we show a comparison between the two codes 
reporting the power spectrum of the central density oscilla- 
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FIG. 18: Central mass-density, normalized to the initial value, in 
a stable TOV star (A/ = 1.4 Af© and polytropic exponent F = 2) 
evolution at different resolutions. PPM and Marquina were used for 
all runs. 



FIG. 19: Power spectrum of the central mass-density evolution of an 
M = 1.4 A/©, F = 2 stable TOV star performed with 128^ grid 
points. The units of the vertical axis are arbitrary. 



tions computed with the Whisky code and the corresponding 
frequencies as obtained with perturbative techniques and with 
the 2D code. Clearly the agreement is very good with an er- 
ror below 1 % in the fundamental frequency. The fact that the 
frequencies computed with the code coincide with the physi- 
cal eigenfrequencies calculated through perturbative analysis 
allows us to study with our code the physical properties of lin- 
ear normal-modes of oscillation even if such oscillations are 
generated numerically. 

The last test performed in the linear regime consists of the 
evolution of stationary solutions of rapidly rotating stars, with 
angular velocity up to 95% of the allowed mass-shedding 
limit for uniformly rotating stars. In this case, a number of 
small improvements on the boundary and gauge conditions 
have allowed us to extend considerably the timescale of our 
evolutions of stable rapidly rotating stars, which can now be 
evolved for about 10 ms, a timescale which is 3 times larger 
than the one previously reported in 1 14]. In analogy with the 
non-rotating case, the truncation error triggers quasi-radial os- 
cillations in the star. Such pulsations converge to zero with 
increasing resolution. Determining the frequency spectrum 
of fully relativistic and rapidly rotating stars is an important 
achievement, allowing the investigation of a parameter space 
which is astrophysically relevant but too difficult to treat with 
current perturbative techniques. More details on this investi- 
gation will be presented in a separate paper 

Finally, we have considered tests of the non-linear dynam- 
ics of isolated spherical relativistic stars. To this purpose we 
have constructed TOV solutions that are placed on the un- 
stable branch of the equilibrium configurations (see inset of 
Fig. I20> . The truncation error in the initial data for a TOV 
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FIG. 20: Normalized central mass-density evolution of an M — 
1.4 Mq, r = 2 unstable TOV star performed with 96"^ grid points. 



solution is sufficient to move the model to a different configu- 
ration and in Whi sky this leads to a rapid migration toward a 
stable configuration of equal rest mass but smaller central den- 
sity. Such a violent expansion produces large amplitude radial 
oscillations in the star that are either at constant amplitude, 
if the polytropic EOS ( I2.12> is used, or are damped through 
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shock heating, if the ideal fluid EOS ( I2.14> is used and the 
equation for t is evolved in time. A summary of this dynam- 
ics is presented in Fig.|20| which shows the time series of the 
normalized central density for a TOV solution. Note that the 
asymptotic central density tends to a value corresponding to a 
rest mass slightly smaller than the initial one (straight dotted 
line). This is the energy loss due to the internal dissipation. 
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